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1.   INTRODUCTION 

Angstroem  (3,  4)  reported  first  measurements  of  thermal  diffusivity 
of  metals  by  the  temperature  wave  method .more  than  100  years  ago.  Since 
that  time  numerous  investigators  have  examined  the  method  as  the  measuring 
equipment  improved  (21,'  30,  29).  Generally,  the  accuracy  of  the  results 
obtained  was  comparable  to  stationary  measurements,  yet  the  method  never 
achieved  widespread  use  since  it  is  more  time  consuming,  both  in  acquisi- 
tion and  interpretation  of  data.   The  method  has  an  advantage  for  the 
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measurement  of  the  thermal  diffusivity  (dimensions:   length  /time)  since 

only  distance  and  time  measurements  are  necessary,  compared  to  the  conduc- 
tivity measurements  where  heat  flow  also  has  to  be  determined. 

In  the  temperature  wave  method  the  measurements  are  made  when  periodic 
temperature  variations  are  achieved  in  the  specimen.   Related  to  this 
method  are  the  transient  methods,  where  the  temperature  development  is 
followed  after  a  pulse  or  step  input  is  applied  to  the  specimen.   A 
transient  method  was  developed  recently  by  several  investigators  (11,  13, 
26)  for  measurements  of  the  thermal  diffusivity  of  small  samples.   The 
temperature  history  of  the  back  side  of  a  thin  sample  is  followed  after 
a  pulse  of  intensive  radiation  is  delivered  to  the  front  side  from  a 
laser  beam.   The  method  is  fast  and  applicable  to  a  wide  range  of  mat- 
erials and  temperature  levels  and  has  been  successfully  used  in  practical 
measurements  (22,  25). 

The  present  investigation  was  in  part  stimulated  by  the  intensive 
investigation  of  neutron  diffusion  waves  (31) ,  which  are  governed  in 
the  diffusion  approximation  by  the  same  equations..  Some  similarities 


exist  although  the  differences  are  probably  greater  than  the  similarities. 
One  of  the  interesting  similarities  is  that  of  the  concept  of  reflection 
and  transmission  of  diffusion  waves.  Since  temperature  waves  (or  heat 
waves)  fall  in  this  category,  it  is  thought  that  the  contact  resistance 
for  the  heat  flow  across  an  interface  might  be  better  understood  with 
the  application  of  this  analytical  technique.  As  no  literature  has  been 
uncovered  which  reports  on  the  use  of  this  technique  for  this  purpose, 
it  may  be  of  great  value.   The  success  of  this  investigation  may  have 
great  practical  application  for  the  measurement  of  the  thermal  contact 
resistance  between  the  cladding  and  the  fuel  in  a  nuclear  fuel  element. 
The  idea  is  to  determine  the  thermal  properties  of  the  assembly  from  the 
response  of  the  cladding  temperature  to  variations  (oscillations)  of  the 
fuel  power  level.   BURDG  et  al.  (9)  performed  some  experiments,  but 
concluded  that  the  sensitivity  of  the  method  is  not  satisfactory,  in  part 
due  to  the  unrefined  data  analysis. 

The  uncertainty  of  the  contact  resistance  in  nuclear  fuel  elements 
for  a  while  spurred  the  research  in  this  field.   Since  nuclear  fuel 
elements  are  currently  designed  by  experience,  good  predictions  are  still 
lacking  for  the  contact  thermal  resistance.   Presently  the  most  interest 
for  this  problem  is  centered  in  the  area  of  space  technology. 

The  additional  resistance  to  the  heat  flow  at  the  interface  of  two 
solids  arises  primarily  from  the  reduced  heat  flow  cross  section  area, 
as  the  real  surfaces  contact  only  on  a  fraction  of  the  apparent  contact 
area.   This  is  due  to  the  waviness  and  roughness  of  the  surfaces.   The 
heat  flow  is  split  at  the  interface  and  constricted  to  individual  points 


of  real  contact.   The  problem  of  predicting  the  contact  resistance  is 
usually  subdivided  into  three  steps:   determination  of  the  constriction 
resistance  for  an  elemental  channel  (region  of  the  solids  associated 
with  a  single  contact  point),  determination  of  the  number  of  contact 
points,  their  dimensions  and  distribution,  and  third,  superposition 
of  the  elemental  channel  solutions.   Early  investigators,  e.g.  Holm 
(13) ,  assumed  rather  idealized  conditions  of  uniformly  spaced  contact 
points.   Later  much  attention  was  given  to  the  determination  of  the 
number  and  distribution  of  the  real  contacts  (10,  32).   Here  important 
influences  are  pressure,  roughness,  waviness,  indentation  hardness  and 
the  elastic  properties  of  the  solids  in  contact.  Most  investigators 
that  were  attempting  to  correlate  the  experiments  with  the  predictions 
did  so  for  nominally  flat  surfaces,  the  specimens  were  prepared  by 
roughening  surfaces  that  were  previously  polished  to  optical  flatness. 
In  machined  joints  the  waviness  is  however  often  more  important  than  the 
roughness.   The  third  task,  superposition  of  the  solutions  for  the 
elemental  channels  of  various  sizes  and  nonuniform  distribution  has  also 
not  yet  been  solved  satisfactorily.   Recently  some  progress  has  been 
made  in  this  area  (24,  32). 

An  interesting  anomaly  in  the  contact  resistance  was  noted,  to  the 
elucidation  of  which  the  temperature  wave  method  may  prove  valuable. 
It  was  noticed  that  in  an  aluminum-steel  pair  the  resistance  is  higher 
for  the  heat  flow  from  steel  to  aluminum  than  in  the  reverse  direction 
(6,  28).   The  difference  amounted  to  approximately  20%  in  air  and  100% 
in  vacuum  being  approximately  the  same  in  absolute  value,  in  the  measure- 


ments  reported  by  Rogers  (28).   Postulated  causes  for  this  effect  were 
attributed  to  the  influence  of  local  strains  (32)  and  to  different 
electric  potentials  of  the  two  metals  that  impede  the  movement  of 
electrons  in  one  direction  (28) .  None  of  the  predicting  theories  includes 
a  satisfactory  explanation  for  this  directional  effect.   With  the  temper- 
ature wave  method,  a  new  mode  of  heat  flow  is  available.   Periodic  heat 
flow  is  applied  and,  more  important,  the  temperature  differences  can  be 
quite  small  compared  to  the  steady  state  measurements. 

Very  few  references  exist  in  the  literature  concerning  the  tempera- 
ture wave  method  of  measuring  the  thermal  resistance  of  contacts.   Green 
(16)  investigated  the  thermal  impedance  of  the  surface  layers  but  his 
work  is  only  partially  relevant  to  this  investigation.   Alzofon  (2) 
theoretically  investigated  the  optical  properties  of  temperature  waves, 
but  with  the  intention  of  constructing  the  steady  state  solutions. 
Heasley  (17)  and  others  (12,  18)  treated  the  transient  contact  resistance 
for  the  evaluation  of  the  heat  transfer  between  bodies  at  different 
temperatures  brought  into  contact.   The  problem  in  this  later  case  is 
symetric  while  when  a  temperature  wave  approaches  a  discontinuity  the 
situation  in  the  two  solids  is  not  identical. 

Based  on  the  reported  work  in  the  areas  briefly  outlined  above  the 
purpose  of  this  investigation  was  set  at  designing  and  building  the 
equipment  for  the  temperature  wave  measurements  and  establishing  the 
method  of  analyzing  the  data.   Both  results  can  be  used  in  the  further 
investigation  of  the  contact  thermal  resistance  and  the  application  of 
the  temperature  wave  method  in  other  research. 


2.  THEORY 

In  this  investigation  the  temperature  waves  were  generated  in  a 
metallic  rod  or  in  an  assembly  of  two  contacting  rods  by  alternating 
the  heat  input  at  one  end  and  keeping  the  other  end  at  a  constant 
temperature.  The  rods  were  of  substantial  diameter  and  well  insulated. 
In  the  case  of  two  contacting  specimens  the  interface  was  perpendicular 
to  the  axis  of  the  specimens.   Thus  the  heat  conduction  problems  are 
of  a  one  dimensional  transient  nature.   The  analytical  solutions  of 
these  problems  are  investigated  in  this  section. 

For  the  case  of  two  contacting  solids  a  constant  heat  transfer 
coefficient  is  considered  to  be  applicable.   It  may  be  questionable 
whether  some  variable  heat  transfer  coefficient  would  be  more  appro- 
priate because  the  heat  transfer  is  periodic.   Also  it  may  be  doubtful 
that  the  heat  transfer  coefficient  so  measured  corresponds  to  the 
stationary  measured  quantity.   In  Appendix  A  the  transient  response  of 
an  elemental  heat  channel  is  investigated  and  it  is  shown  that,  by 
considering  the  dimensions  of  the  real  contacting  points  and  their 
spacing  in  the  contacts  of  flat  metallic  surfaces,  the  time  scale  on 
which  the  transient  effects  are  important  is  of  the  order  of  .1  sec. 
The  periods  considered  for  the  experiments  are  of  the  order  of  100 
seconds  and  thus  the  transient  effects  certainly  can  be  neglected. 

Whether  the  heat  transfer  coefficient  measured  by  the  periodic 
heat  flow  method  is  identical  to  the  one  obtained  by  a  stationary 
method  cannot  be  resolved  at  this  time.   A  possible  reason  for  a 
systematic  difference  could  be  a  directional  heat  flow  effect.   For 


steel-aluminum  contacts,  differences  of  up  to  100%  have  been  observed 
between  measurements  with  different  directions  of  heat  flow  (16,  17). 
As  the  exact  mechanism  of  this  effect  is  not  known,  a  similar  discre- 
pancy for  the  periodic  heat  flow  case  cannot  be  ruled  out  but  should 
be  investigated  experimentally. 

2.1   Temperature  Waves  in  a  Slab 

2.1.1  Stationary  Eoundary  at  Constant  Temperature 

The  one  dimensional  heat  conduction  problem  is  solved  for  a  slab 
with  periodic  temperature  on  one  boundary  and  constant  temperature  on 
the  other  boundary.  The  governing  equation  is 


32T(x,t)  =  3T(x,t) 

„  2     at 
ax 


(i) 


where  a  =  —. 

PCP 


The  initial  condition  is 


T(x,0)  =  0  (2) 

and  the  boundary  conditions  are 

T(L,t)  '=  A  sin  ut  t  >  0 

T(0,t)  =  0  t  >  0  (3) 

To  obtain  the  solution  of  this  problem  the  Laplace  transform  technique 


is  used.   Transforming  the  variable  t  into  the  Laplace  variable  s 
the  governing  equation  is 


2- 

a   d  T(^'S)  =  sT(x,s)  -  T(x,0) 
dx 


applying  the  initial  conditions  this  reduces  to 


(4) 


2- 

d  T(x,s)   s  -,   .  ,,  . 
V  =  -  T(x,s)  (4a) 

dx2    a 


This  equation  has  the  general  solution 


ICx.b)  «  CjfeK"   +  C2(s)e  ,uc  (5) 

The  transformed  boundary  conditions  are 

T(0,s)  =  0 

T(L,s)  =  -j^-j  (6) 

u  +  s 


From  these  boundary  conditions  the  constants  C.  and  C„  are  evaluated: 


C2  "  "  Cl 


Aw 


Cl   sinh  (/s  L)(co2  +  s2)  (7) 


,    sinh  f—  x 

Kx.s)  •  -jj 2  77~~  () 

d)  +  s  sinh  /L  L 


!  a 


The  inverse  transform  is  found  by  the  theorem  of  residues,  as  the  only 
singularities  are  poles.   Furthermore  all  poles  are  simple.   The  first 
part  of  the  denominator  has  two  zeroes,  s  =  +  iw.   It  can  be  shown. 
that  the  only  zeros  of  the  function  sinh  |—  L  are  at 


—  L  =  m  tt 
la 


and  thus 


2  2 
^a  n-  1,  2,  ...  (9) 

L 


Omitting  the  algebraic  steps  in  the  evaluation  of  the  residues,  the 
final  result  is 


Residue  at  (iw)  +  Residue  at  (-ito)  ■ : — s—r: n^r   [sin  cot 

cosh  2AL  -  cos  2AL 


(sinh  A  x  cos  A  x  sinh  A  L  cosh  A  L  +  cosh  A  x  sin  A  x  cosh  A  L  sin  A  L) 


4-  cos  tot  (-sinh  A  x  cos  A  x  cosh  A  L  sin  A  L  +  cosh  A  x  sin  A  x 


sinh  A  L  cos  A  L) ] 


where   A  ■ 


fe  <10> 


2 
The  residue  at  any  pole,   s  = x—  a,   is 


L2 


2 

Residue  at  (-  - — « — ) 
L 


2  2 

+1  an  tt 

2Aiu(-l)  ' .   mix  "  L2  ,,,% 

, ,      —IT"  T e  (11) 

(A2  +  a2  HJL_  ) 
IT 


These  terms  decay  with  time  and  represent  the  transient  part  of  the 
solution,  while  the  terms  corresponding  to  the  poles  +  iu  represent  the 
asymptotic  part  of  the  solution.   The  transient  solution  is  not  partic- 
ularly interesting  as  the  measurements  are  performed  after  steady  state 
is  achieved.   The  magnitude  of  the  transient  part  is  briefly  evaluated 
as  follows. 

2lT  Jlrt 

Introducing  t  =  —  (period  of  oscillation)  and  L,  ,  «  1/A  =  I  — 
o    a  1/e         I  ui 

(attenuation  distance)  the  time  necessary  for  the  slowest  decaying 
term  to  be  reduced  by  a  certain  factor  can  be  evaluated.   In  Table  I 
the  number  of  periods  for  a  100  fold  reduction  of  the  slowest  decaying 
term  in  slabs  of  different  thicknesses  are  presented. 

Table  I 

Evaluation  of  the  transient  solution. 

Dimensionless  thickness   Reduction  of  the  first      No.  of  periods  for  a 
of  the  slab,  I/I-,/       transient  term  per  period    100  fold  reduction 

4  .135  2.2 

6  .422  5.2 

8  .618  9.5 

10  .734  14.5 
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2.1.2  Asymptotic  Solution  for  a  Thick  Slab. 


Two  useful  approximations  can  be  derived  from  the  asymptotic 
solution  Eq.  (10)  for  thick  slabs  (L/L  .  large).   For  AL  >>  1  the 

approximation 

coshAL  "-'  sinhAL  ■  -r-  e 


can  be  used  and  cos2AL,  which  is  always  less  than  unity  can  be  neglected 
in  the  denominator  of  Eq.  (10).   Substituting  the  exponential  expressions 
for  the  hyperbolic  functions  and  implementing  the  addition  theorems 
for  the  harmonic  functions  the  result 


^(X.t)  =  -rj  [eXX  sin(tot-A(L-x)j  -  e  Ax  sin(oit  -  A(xfL))]   (12) 

e 


is  obtained. 

This  expression  shall  be  called  the  "one  reflection  approximation". 
If  furthermore  x  is  large,  i.e.  only  the  regions  far  from  the 
stationary  boundary  are  observed,  the  result  is 


QT(x,t)  =  Ae"X(L  X)sin(u,t  -  A(L-x))  .  (13) 


This  expression  shall  be  called  the  "no  reflection  approximation". 

The  one  reflection  approximation,  Eq.  (12),  can  be  interpreted  as 
the  sum  of  two  waves,  one  approaching  the  boundary  with  the  constant 
temperature  and  the  other  wave  leaving  this  boundary  with  reversed  sign. 


11. 


Each  wave  has  the  form 

Ce-^  8tn(ast-X£)  (14)' 

where  £  is  the  linear  dimension  measured  in  the  direction  of  the  wave 
travel. 

The  amplitudes  and  the  phase  lags  were  calculated  from  the  three 
expressions,  the  exact  asymptotic  solution,  Eq.  (10),  the  no-  and  one- 
reflection  approximations,  Eqs.  (12)  and  (13);  for  slabs  of  various 
thicknesses.   The  results  are  presented  in  Figures  1  to  3.   From  the 
numerical  results  the  following  can  be  concluded:   For  slabs  thicker 
than  3  Lj.  the  no-reflection  approximation  is  adequate,  except  closer 
than  1.5  L-we  from  the  stationary  boundary  where  the  one  reflection 

approximation  is  to  be  used.   For  slabs  thinner  than  1.5  L,  ,   the  exact 

1/e 

solution  is  to  be  used  instead  of  the  reflection  approximations. 
2.1.3.   Stationary  Boundary  Insulated 

The  governing  equation,  initial  condition  and  the  periodic  boundary 
condition  are  same  as  before,  Eqs.  (1),  (2)  and  (3),  but  on  the  stationary 
boundary  the  following  condition  is  applied: 


I  ».*)  "  °  (15) 
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Exact  solution 


One  reflection  approx. 

_  No  reflection  approx. 
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.  Figure  1.   Amplitude  and  phase  of  the  temperature  wave, 

thickness  of  the  slab  »  1  L,  ,  . 

1/e 


13 


Exact  solution  and  one 
reflection  approx. 


No  reflection  approx. 


Oscillatory 
Boundary 


Ax 


Stationary 
Boundary 


Figure  2.  Amplitude  and  phase  of  the  temperature  wave,  thickness 

of  the  slab  =  3L_  ,    . 
1/e 
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Figure  3.  Amplitude  and  phase  of  the  temperature  wave, 

thickness  of  the  slab  =  5  L,  , 
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Due  to  this  boundary  condition  the  relation  of  the  constants  in  the 
Laplace  transformed  solution  is 


^(s)  =  C2(s)  (16) 


Further  the  solution  is  obtained  in  a  similar  way  as  before,  and  the 
result  is 

.     cosh  [—  x 

*<*■">  "-rr-T  — ■ >-  <17> 

10  +  s   cosh  I—  L 


The  inverse  transform  can  again  be  obtained  by  evaluation  of  the  resi- 
dues.  The  first  part  of  the  transform  has  two  simple  poles,  s1,„=+iu 

and  it  can  be  shown  that  the  function  cosh  f—  L  has  zeroes  at 

la 


(2n-l)2n2 

2  2 
2  L 


(18) 


i.e.  when 


/s  ,    .  ,2n-l.  , 

|-L  =  i(— 2~>  *  n  =  1,  2,  .... 

These  poles  are  on  the  negative  real  axis  and  thus  the  contribution 
of  their  residuals  is  of  transient  nature.  The  residues  at  the  two 
poles  on  the  imaginary  axis  form  the  asymptotic  solution; 
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Residue  at  (iw)  +  Residue  at  (-iu)  = 

=  7 ,  „.  — ; ^TTV  [cosut(-coshXxcosXXsinhXL  sinXL  + 

(cosh/AL  +  coszXL) 

+  sinhX x sinX x  coshXL  cosXL)  +  sinwt (coshX x  cosXx  coshXL  cosXL  + 

+  sinhX  x  sinX  x  sinhXL  sinXL) ]  (19) 

For  thick  slabs  using  the  same  assumptions  as  in  section  2.1.2,  the 
result  becomes 


xT(x,t)  =  Ae  XL[eXXsin(iot-Xa-x))  +  e~Xx  sin(iot-X  (L+x))  ]        (20) 


This  result  can  be  interpreted  as  consisting  of  two  parts,  a  first  term 
representing  a  wave  approaching  the  insulated  boundary  at  x  =  0  and  a 
second  term  representing  a  wave  leaving  this  boundary.   At  the  boundary 
the  two  waves  are  in  phase  and  reinforce  each  other,  while  in  the  case 
of  the  constant  boundary  temperature  they  were  of  opposite  phase. 

2.2  Temperature  Wave  at  a  Discontinuity 

The  behavior  of  a  plane  temperature  wave  at  a  contact  of  two  solids, 
the  contacting  surface  being  perpendicular  to  the  direction  of  the  wave 
travel,  shall  be  investigated.   The  heat  transfer  characteristics  of 
the  discontinuity  is  assumed  to  be 

q"  -  hAT  •  (21) 
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where  q"  is  the  heat  flow  per  unit  area,  h  is  the  heat  transfer  coef- 
ficient and  AT  is  the  temperature  drop  across  the  discontinuity.   The 
coordinate  x  is  placed  parallel  to  the  direction  of  the  approaching 
temperature  wave  (perpendicular  to  the  discontinuity)  and  the  origin  is 
placed  at  the  discontinuity.   In  analogy  with  the  case  of  the  insulated 
boundary  or  the  stationary  boundary  the  solution  is  expected  to  consist 
of  a  reflected  wave  with  some  phase  change  but  with  the  same  period  and 
a  transmitted  wave,  also  with  some  phase  change.  The  forms  assumed  are, 
the  amplitudes  being  normalized  so  that  the  approaching  wave  has  ampli- 
tude 1  at  discontinuity: 
Approaching  wave 


—A  x 
e  1  cos(ut  -  A.x)  (22) 


Reflected  wave 


A  x 
Ar  e   1     cos(ut  +  6r   +  A.x)  (23) 


Transmitted  wave 


—A  x 
At  e  2  COB C«t  +  &t   -  A2x)  (24) 


Taking  the  amplitude  of  the  approaching  wave  as  unity  at  x  =  0,  the  time 
dependent  temperatures  in  solids  1  and  2  are 
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^(x.t)   =  e  Xlx  cos  (ut-A  ax)  +  A  eAlX  cosfot  +&     +  X.x)      x  <   0        (25) 


-A„x 


T2(x,t)   »  Afce     2      cos(ut  +5t-A2x)  x  >   0  (26) 

Substituting   these   expressions   into  the  equations   of  heat  balance  at   the 
discontinuity 


9T 
"  kl  23T  (0,t)   =  5"(0,t)    =  h(T1(0,t)    -  T2(0,t))  (27) 

3T 
"  k2   3x~~  (0'°   =  q''^0'^   "  h(T1(0,t)    -  T2(0,t))  (28) 


yields 

-  k^A-^A    (cos(iut+(5    )   -  sin(uit+S    ))   -    (costot-sinut)]   = 

■  h[cosut  +  A  cos(ut+5    )    -  A  cos(wt+i5    )]  (29) 

and. 

k2*2At[cos((i)t+6    )    -  sln(iat+S    )]    = 

=  hlcosut  +  A  cos(ut+S  )  -  A  cos  (uit+5  )]  (30) 

The  material  constants  and  the  heat  transfer  coefficient  can  be  grouped 
in  dimensionless  quantities 

Yx  -  h/k1A1        y2  *  h/k2A2  (31) 
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This  quantity  relates  the  heat  transfer  distance  k/h,  i.e.  the  thickness 

of  the  solid  causing  the  same  temperature  drop  as  the  discontinuity, 

/2 — 

to  the  1/e  length  of  the  temperature  waves,  L  ,  =  1/A  =  —  . 

1/e  to 

After  using  the  summation  theorems   for   the   trigonometric  functions 
and  observing  that  the  identity  has   to  be  conserved  at  all  times,   four 
equations  are  obtained  from  Eq.    (29)    and    (30) 


-(1+Y,)  A  cos<5     +  A  sinS     +  Y,A  cos<$  »  -Q-Y, ) 

lrr  rr  ltt  1 

Arcos<5r  +   (l-Y1)ArsinSr  -  Y^sin^  =   1 


-Y^,-003^  +  U+Y2)Atcos6t  -       A  sin« 


2 


Y2Arsin6r  +  Atcos6t+(l+Y2)  Atsin6t   =  0 

(32) 

From  this  system  of  linear  equations  for  A  cosS  ,  A  sinS  ,  A  cos6  and 

r    r   r    r   t    t 

A  sinS  these  quantities  can  be  calculated  for  any  values  y     and  y 
and  from  them  the  amplitude  and  phase  of  the  temperature  oscillation 
at  any  point  in  both  solids  can  be  obtained. 

The  asymptotic  solution  for  y     and  y     approaching  infinity  cannot 
be  obtained  from  the  system  Eq.  (32),  as  the  problem  formulation  Eqs. 
(27),  (28)  is  not  correct  in  this  case.   The  following  conditions  have 
to  be  used  at  the  interface 

TjCO.t)  -  T2(0,t)  (33) 
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3T  3T 

k.  T~   (0,t)  =  k,  -zf-   (0,t) 


1  3x 


2  3x 


(34) 


Introducing  the  assumed  wave  forms,  Eqs.  (22),  (23),  (24),  performing 
the  algebraic  manipulations  and  solving  the  resulting  linear  system  of 
equations  the  following  is  obtained 


At  cos6t  -  §L 


A  cos  5 


r-i 
r+i 


A  sinS 


A  sin5 


(35) 


where  r  is  defined  as 


k2x2 


/kicPipi  V; 


k2Cp2P2   Yl 


(36) 


The  quantities  A ,  5  ,  A  ,  6   as  functions  of  y      are  plotted  for 
various  ratios  -/2^l   in  F18ures  4  t0  ?•   The  asymptotic  values  are 
also  indicated.   The  value  of  the  ratio  y  /y  ,  i.e.  r,  for  metallic 
pairs  ranges  between  4  and  1/4.   Of  all  metals  the  best  conductor- 
capacitor  is  copper,  with  the  kp c  value  of  1,390  J2/m*c2sec,  while 
poor  conductors  and  capacitors  like  stainless  steel  have  values  around 
80  for  this  parameter. 
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The  graphs  of  the  phase  angles  and  the  reflected  and  transmitted 
fractions  reveal  some  interesting  properties  (Figures  5  and  7).   The 
reflected  wave  always  preceeds  the  incoming  wave  (positive  angle  6  ) 
and  the  angle  asymptotically  approaches  three  values,  0,1/4  and  H, 
depending  on  the  tJy1   ratio.   For  T   less  than  one,  i.e.  when  the 
second  material  is  a  better  conductor- capacitor,  the  reflected  wave 
has  opposite  phase  to  that  of  the  incoming  wave.   The  same  relation 
was  obtained  in  the  case  of  the  stationary  boundary  temperature  case, 
only  now  the  amplitude  is  less  than  that  of  the  incoming  wave  and 
depends  on  the  r  value.   For  r  larger  than  one  the  phase  approaches 
zero  for  large  Yj  and  the  situation  is  similar  to  the  insulated  boundary 
case.  Note  also  that  for  r  larger  than  one  the  transmitted  wave  ampli- 
tude rises  above  unity,  i.e.  is  larger  than  that  of  the  approaching 

wave.  This  does  not  violate  the  boundary  conditions  at  the  interface 

as  the  reflected  wave  is  in  phase  with  the  incoming  wave  and  the  sum 
(vectorial)  is  always  larger  than  the  amplitude  of  the  transmitted 

wave. 

For  the  ideal  contact  (y-«,  h-«°)  only  two  values  of  the  reflected 

wave  phase  angle  are  possible,  0  and  H.   The  curve  f or r  =1  approaches 

Jl/4,  but  at  infinity  the  amplitude  is  zero. 

In  Fig.  8  the  phase  difference  between  the  sum  of  the  approaching 

and  the  reflected  wave  and  the  transmitted  wave  at  the  interface  is 

plotted. 

In  Fig.  9  to  Fig.  12  the  phase  lag  profiles  in  the  incoming  solid 

are  presented.   Note  that  here  the  phase  lag  is  plotted,  which  is 

defined  with  the  opposite  sign  as  the  phase  difference. 
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3.  EXPERIMENT  DESIGN 

The  purpose  of  the  experimental  work  was  to  review  the  temperature 
wave  method  of  thermal  diffusivity  measurements  and  to  investigate  the 
feasibility  of  contact  resistance  measurements  by  this  method.   Thus 
there  were  two  groups  of  preceeding  experiments  to  be  considered,  the 
temperature  wave  measurements  and  the  contact  conductance  measurements. 
To  accomplish  the  combined  task  the  present  experiment  differs  in  some 
respects  from  both  groups  of  prototypes. 

Angstroem  (3)  and  most  later  investigators  of  the  temperature  wave 
method  (21,  29,  30)  did  not  attempt  to  restrict  the  surface  heat  losses 
from  the  specimens,  which  were  in  some  cases  in  the  form  of  wires.   They 
rather  accounted  for  these  losses  in  the  analysis.   In  this  case  the 
temperature  oscillations  have  to  be  measured  at  least  at  three  points. 
For  the  contact  resistance  measurement  the  specimens  have  to  be  of 
substantial  cross  section  and  the  heat  losses  would  make  the  temperature 
field  two  dimensional  due  to  the  radial  flow.   To  avoid  any  such  further 
complication,  as  the  contact  between  the  two  specimens  already  introduces 
some,  it  was  considered  preferable  in  the  present  experiment  to  eliminate 
these  surface  heat  losses  if  possible. 

In  the  contact  resistance  experiments,  two  specimens  are  usually 
sandwiched  between  calibrated  heat  flow  metering  sections,  possibly 
other  layers  of  different  materials,  with  the  heater  and  cooler  at 
opposite  sides.   In  such  arrangements  the  temperature  waves  would  be 
greatly  attenuated  at  each  interface.   For  this  reason  the  specimens 
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should  be  in  direct  contact  with  the  heater  and  the  cooler..   The  surface 
heat  losses  were  eliminated  in  the  contact  resistance  and  conductivity 
measurements  by  placing  a  guard  sleeve  around  the  specimen  with  the  same 
temperature  gradient  as  in  the  specimen.   Due  to  the  temperature  oscilla- 
tion, only  the  mean  temperature  can  be  matched  in  this  way.   To  avoid 
the  heat  losses  due  to  the  deviations  from  this  mean,  insulation  is 
required. 

The  basic  differences  of  the  present  design  from  the  two  groups  of 
prototype  experiments  are  that  the  specimens  are  of  substantial  diameter 
and  are  insulated  and  in  good  contact  with  the  heater  and  cooler.   The 
length  of  the  specimens  was  chosen  primarily  with  regard  to  the  precision 
of  the  thermocouple  positioning  and  the  thermocouple  dimensions.   The 
temperatures  were  measured  on  the  axis  of  the  specimens  and  the  thermo- 
couples placed  in  holes  drilled  into  the  specimens.   The  position  of  the 
thermocouple  in  the  hole  is  not  well  determined,  so  the  diameter  of  the 
hole  must  be  considered  as  the  uncertainty  interval  for  the  thermocouple 
location.   With  the  hole  diameters  used,  .03"  or  .05",  and  the  thermo- 
couple separation  distance  of  2",  the  uncertainty  is  3%  or  5%,  respectively, 
and  it  was  considered  that  a  shorter  distance  should  not  be  used.   The 
length  of  the  contact  resistance  specimens  was  7  1/2"  and  the  length  of 
the  thermal  diffusivity  measurement  specimens  was  15".   The  periods  of 
oscillation  to  be  used  were  determined  from  the  distance  between  the 
temperature  measuring  points,  so  that  the  attenuation  would  not  be  exces- 
sive.  If  the  amplitude  of  the  oscillation  at  the  first  thermocouple  is 
of  the  order  of  a  few  degrees  (higher  amplitude  would  be  undesirable) , 
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and  the  estimated  limit  of  detection  is  of  the  order  of  .01  C,  it  was 
determined  that  the  distance  between  the  extreme  positions  should  not 
exceed  5  attenuation  lengths,  i.e.  the  attenuation  length  should  not 
be  shorter  than  1/5  of  the  longest  distance. 

In  the  first  design,  the  heater,  which  generates  the  temperature 
oscillations  of  one  end  of  the  specimen,  was  equipped  with  a  copper 
plate  at  the  bottom  which  was  placed  in  contact  with  the  end  of  the 
specimen  rod.   Such  arrangement  produced  very  small  oscillations  of 
temperature  in  the  specimen  due  to  relatively  large  thermal  capacity 
of  the  heater  and  poor  contact  with  the  specimen.   The  design  was 
subsequently  changed  in  favor  of  a  band  heater  tightened  around  the 
end  of  the  specimen.   In  this  case  doubts  arise  concerning  the  uniformity 
of  the  temperature  across  the  specimen  bar.   A  numerical  analysis  of 
the  temperature  profile  was  performed  and  in  Fig.  13,  the  lines  of 
equal  temperature  are  plotted  for  a  rod  with  constant  temperature  on 
part  of  its  circumference  and  uniform  heat  flow  at  a  distance  from  the 
end.   It  is  evident  that  the  heat  flow  becomes  completely  one  dimensional 
less  than  one  diameter  from  the  limits  of  the  uniform  temperature  of  the 
perimeter.   It  is  assumed  that  the  same  is  valid  for  the  periodic  heat 
flow  and  a  conservative  margin  of  two  rod  diameters  from  the  heater  was 
used  for  the  positioning  of  the  first  thermocouple.   The  rod  diameter  of 
one  inch  was  chosen  arbitrarily. 

This  arrangement  produced  temperature  oscillations  of  sufficient 
amplitude,  hut  a  high  temperature  gradient  was  needed  to  deliver  the 
heat  to  the  cooler  at  the  other  end  of  the  specimen.   Thus  a  potentially 
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Figure  13.   Temperature  distribution  in  a  cylindrical  rod  with  constant 
circumferential  surface  temperature  at  one  end. 
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important  advantage  of  the  temperature  wave  method,  that  it  does  not 
intrinsically  require  large  temperature  gradients,  was  lost.   Subse- 
quently a  cooler  was  installed  at  this  end,  and  the  flow  of  coolant  was 
oscillated  at  the  same  frequency  as  the  heater.   However,  the  solenoid 
valve  used  caused  unacceptable  disturbances  in  the  thermocouple  voltage 
output  due  to  magnetic  induction,  so  the  valve  was  removed  and  the 
coolant  flow  was  maintained  constant.   Contrary  to  the  expectations, 
the  temperature  oscillation  amplitude  did  not  decrease  significantly 
when  the  cooler  was  operated  continuously,  and  the  constant  temperature 
level  decreased.  This  small  effect  of  coolant  flow  can  be  attributed 
to  poor  heat  transfer  from  the  specimen  to  the  coolant  as  the  coolant 
velocity  was  quite  low. 

The  theory  of  the  temperature  waves  was  developed  for  sinusoidal 
temperature  variation  on  the  boundary.  To  approximate  this  condition 
a  mechanism  was  constructed  to  vary  the  voltage  on  the  heater  harmoni- 
cally.  It  was  soon  percieved  that  this  is  not  necessary,  as  any  higher 
harmonic  components  are  very  weakly  excited  in  the  specimen  and  die  out 
rapidly  along  its  length.  A  cam  operated  on-off  control  of  the  heater 
was  installed.   Kith  this  arrangement,  the  first  higher  harmonic  was 
detected  only  in  the  first  thermocouple  output  and  had  a  barely  signi- 
ficant amplitude. 

The  largest  body  of  thermal  contact  conductance  data  available  is 
for  contacts  in  vacuum.   The  specimen  and  heater  housing  in  this  investi- 
gation- was  designed  for  vacuum,  but  due  to  undetected  faults  the  required 
vacuum  was  not  achieved. 


36 


Due  to  the  inherent  attenuation  characteristics  of  the  temperature 
waves  it  was  anticipated  that  the  signal  to  noise  ratio  from  the  more 
distant  thermocouples  would  require  a  refined  data  analysis  technique,  , 
considering  the  accuracy  required.   The  correlation  technique  offers 
a  possibility  of  analyzing  noisy  signals.   An  outline  of  the  technique 
is  given  in  Appendix  B. 

The  date  acquisition  and  analysis  equipment  is  centered  around 
the  correlation  computer.   The  thermocouple  outputs  were  first  recorded 
on  a  magnetic  tape  recorder  and  analyzed  later.   A  few  initial  measure- 
ments were  performed  by  recording  the  temperature  oscillations  on  a 
strip  chart  recorder  and  analyzing  the  records  manually,  but  this  proved 
highly  unpractical.   One  problem  of  the  temperature  measurements  was 
associated  with  the  elimination  of  the  electronic  noise  that  was  picked 
up  by  the  thermocouple  leads.   The  noise  was  mostly  of  high  frequency 
and  was  eliminated  by  passing  the  signal  through  a  low-pass  filter  before 
recording.   Later  it  was  found  that  in  the  tape  recorder  itself  some  high 
frequency  noise  is  generated  and  another  low-pass  filter  was  used  between 
the  tape  recorder  and  the  correlator.   The  analysis  of  these  circuits  is 
given  in  Appendix  D. 

The  description  of  the  present  equipment  is  given  in  the  following 
section  and  an  evaluation  in  section  8. 
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4.   EXPERIMENTAL  APPARATUS 

The  experimental  equipment  consists  of  the  mechanical  and  thermal 
system  and  of  the  data  acquisition  system.   The  individual  parts  of 
the  equipment  are  described  below  from  the  specimens  to  the  correlation 
computer. 

4. 1  Specimens 

For  the  two  types  of  measurements,  the  thermal  diffusivity  measure- 
ments and  for  the  contact  conductance  measurements,  two  types  of  specimens 
were  used.  The  contact  resistance  specimens  were  used  in  pairs,  the 
total  length  of  the  pair  being  the  same  as  a  single  diffusivity  specimen. 
The  choice  of  the  dimensions  of  the  specimens  is  discussed  in  section  3. 

4.1.1  Thermal  Diffusivity  Specimens 

Solid  rods  of  the  material  under  consideration,  1  inch  in  diameter 
and  15  inches  long,  were  used  for  the  thermal  diffusivity  measurements. 
Five  holes  were  drilled  in  each  specimen,  2  1/2  inches  apart,  the  first 
hole  being  3  inches  from  the  end  that  was  to  be  subjected  to  periodic 
temperature  conditions.   The  holes  were  .05  inch  diameter  and  1/2  inch 
deep,  i.e.  they  reached  the  center  of  the  rod.   Holes  were  drilled  in 
both  ends  of  the  specimen  with  a  #5  countersink  drill  to  a  depth  of 
1/2  inch,  as  required  for  the  cooler  centering.   Thirty  gage  copper 
constantan  thermocouple  wire  was  used  for  the  thermocouples.   The 
junctions  were  formed  by  mercury  pool  immersion  welding  and  cemented 
in  the  holes  with  Sauereisen  cement.   Before  the  installation  of  the 
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thermocouples  the  distance  between  the  thermocouple  holes  was  measured 
precisely. 

A. 1.2  Thermal  Contact  Resistance  Specimens 

The  thermal  contact  resistance  specimens  were  1  inch  O.D.  rods, 
7  1/2  inches  long  (Fig.  14).   Three  holes  were  drilled  along  the  rod 
2  inches  apart,  the  distance  of  the  closest  point  from  the  contacting 
surface  being  1/2  inch.   All  holes  were  .031  inch  diameter  and  1/2  inch 
deep,  the  bottom,  located  at  the  center  of  the  rod.   In  one  end  of  the 
specimen  a  hole  was  drilled  with  a  #5  countersink  drill,  as  required 
for  the  cooler  positioning.   The  other  end  was  ground  to  the  best  possible 
flatness  obtainable  on  a  horizontal  grinding  machine.   After  machining, 
the  contacting  surface  was  carefully  protected.   The  roughness  was 
measured  on  a  Micrometrical  Manufacturing  Co.   Profilometer  (Model  741, 
IE  838)  and  the  hardness  of  the  material  was  tested  on  a  side  surface 
using  the  Rockwell  ball  method. 

Junctions  on  a  36  gage  copper  constantan  thermocouple  wire  were 
made  by  the  mercury  immersion  technique  and  the  thermocouples  were 
cemented  in  respective  holes  with  Sauereisen  cement.   The  electric 
resistance  between  the  two  wires  and  between  the  wires  and  the  specimen 
was  measured  in  order  to  check  the  wire  junction  and  the  junction  contact 
with  the  specimen.  A  photograph  of  the  specimen  is  presented  in  Fig.  15. 

4.2  Mechanical  and  Thermal  Equipment 

The  experimental  vessel  and  all  adjoining  mechanical  and  thermal  - 
equipaient  are  mounted  for  flexibility  of  arrangement  on  a  Unistrut 
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COUNTERSINK    DRILL      9    HOLE 


Figure  14.      Contact   conductance   specimen;    dimensions 
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frame  with  roller  casters.   The  principal  parts  of  the  mechanical  and 
thermal  equipment  are  the  experimental  vessel  with  the  internals ,  the 
loading  mechanism  and  the  heating  and  cooling  control.   The  arrangement 
of  the  mechanical  components  is  shown  in  the  photograph  in  Fig.  16. 

4.2.1  Experimental  Vessel  and  Internals 

The  cross  section  drawing  of  the  experimental  vessel  is  presented 
in  Fig.  17  and  a  photograph  of  internals  in  Fig.  18. 

The  specimen  or  specimen  pair  is  inserted  in  the  lower  cooler,  which 
is  supported  flexibly  on  a  steel  sphere.   The  steel  column  supporting 
the  lower  cooler  is  equipped  with  an  auxiliary  heater  to  prevent  conden- 
sation on  the  cold  surfaces  of  the  lower  section  of  the  experimental 
vessel. 

The  upper  part  of  the  specimen  (or  specimen  pair)  is  equipped  with 
a  band  heater,  rated  at  100  H  at  110  V,  produced  by  Ogden  Manufacturing 
Company.   This  heater  was  controlled  by  a  cam  operated  mlcroswitch  and 
produced  temperature  oscillations  in  the  specimen.  A  force  was  applied 
to  the  top  cooler  by  the  loading  ram,  which  is  equipped  with  strain 
gages  and  is  considered  a  part  of  the  loading  mechanism. 

Both  the  top  and  the  bottom  cooler  tube  connections  are  in  the 
bottom  plate.   The  tubes  coming  out  of  the  bottom  plate  are  copper  with 
Tygon  plastic  connections  for  greater  flexibility.   The  coolers  are 
provided  with  0-ring  seals  to  make  them  water  tight. 

The  heat  losses  from  the  specimens  are  prevented  by  insulation  of 
the  specimens  and  by  the  guard  sleeve,  which  was  maintained  at  the 
same  temperature  as  the  specimens  during  the  measurements .   The  temper- 
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Figure  17.   Experimental  vessel  cross  section 
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Figure  18.   Experimental  vessel  internals. 
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ature  gradient  was  controlled  by  a  heater  at  the  top  end  of  the  guard 
sleeve  and  by  a  cooling  coil  at  the  bottom  end.   Six  thermocouples 
were  installed  in  the  sleeve  and  additional  insulation  was  placed  on 
the  outside.   The  specimen  insulation  is  a  glass  fiber  material  and 
the  guard  sleeve  insulation  a  refractory  felt. 

The  cylindrical  part  of  the  vessel,  made  of  aluminum,  was  bolted 
to  the  bottom  and  top  plates.   All  connecting  surfaces  are  provided 
with  0-ring  grooves,  but  0-rings  were  not  in  place  during  experiments 
as  no  experiment  with  reduced  pressure  was  performed.   The  power  and 
thermocouple  wires  are  conducted  from  the  vessel  through  a  wax  filled 
pot  to  provide  a  gas  tight  seal.   The  loading  mechanism  consists  of 
the  loading  lever  with  the  weight  basket  and  the  load  ram  with  guides. 
The  lower  part  of  the  loading  ram  serves  as  a  load  cell.   At  that 
position  a  hole  is  drilled  in  the  center  of  the  ram  to  allow  larger 
strain.  Two  strain  gages  were  placed  on  the  opposite  sides  of  the  load 
cell  section.  The  strain  gages  were  connected  to  an  universal  amplifier 
with  a  strain  gage  input  box  (Brush  Electronic,  Universal  Amplifier, 
Model  BL-520) .   A  precision  voltmeter  provided  the  output  reading 
(Keithley  Electrometer  600  A,  NE  906).   The  load  cell  was  calibrated 
against  a  load  ring  certified  by  the  National  Bureau  of  Standards.   A 
photograph  of  the  load  cell  portion  of  the  loading  ram  is  presented  in 
Tig..  19. 

4.2.2  Heating  and  Cooling  Control 

The  heating  and  cooling  control  devices  are  located  on  the  control 
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Figure  19.   Load  cell. 
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panel.   Two  autotransformers  for  continuous  power  setting  of  the  heaters 
are  provided,  one  for  the  specimen  heater  and  the  other  for  the  guard 
sleeve  heater.   Both  were  connected  to  unregulated  110  V  mains.  A 
microswitch,  operated  by  a  cam  drive  is  installed  in  the  top  heater 
transformer  power  line.   The  cam  is  driven  by  a  10  rpm  synchronous 
motor  through  a  reduction  gear.   The  reduction  gear  consists  of  a  gear 
pair  and  a  multiratio  gearbox  (Precision  Instrument  Co.,  Dial-A-Ratio 
Gearbox).   The  gear  pair  ratios  used  were  1.3:  1,1:  1  and  1:  1.3  and 
the  available  gearbox  reduction  ratios  were  2.5,  5,  10,  25,  50  and 
100.   Fifteen  different  periods  could  be  obtained  from  .192  to  13  minutes 
with  a  30%  maximum  step  between  two  ratios. 

Cold  tapwater  was  used  as  the  coolant.   The  water  flow  to  the  three 
coolers,  the  top  and  bottom  specimen  coolers  and  the  guard  sleeve  cooler, 
the  top  and  bottom  specimen  coolers  and  the  guard  sleeve  cooler,  was 
monitored  on  three  flow  meters  and  regulated  by  manual  valves.  A 
solenoid  valve  was  provided  for  the  intermittent  operation  of  the  top 
specimen  cooler  but  was  not  used  in  the  experiments. 

4.3  Data  Acquisition  Equipment 

The  specimen  thermocouple  voltage  outputs  were  amplified  and  recorded 
on  a  magnetic  tape  recorder.   Later  the  records  were  analyzed  with  the 
aid  of  a  correlation  computer.   The  correlation  computer  output  represents 
the  raw  data  that  was  further  analyzed  by  the  digital  computer.   For 
steady  temperature  measurements  a  potentiometer  was  used  (Leeds  and 
Northrup,  Portable  Potentiometer). 
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4.3.1  Recording  Equipment 

The  recording  equipment  arrangement  is  shown  in  Fig.  20. 

The  thermocouple  reference  point  was  provided  by  a  thermos  vessel 
filled  with  an  ice-water  slurry.   The  thermocouple  voltage  output  was 
fed  to  the  amplifying  and  filtering  circuits,  composed  of  the  compon- 
ents of  an  analog  computer  (Electronic  Associates  Co.  PACE  TR-10,  NE 
916) .   The  schematic  diagram  and  analysis  of  these  circuits  are  given 
in  Appendix  D. 

The  number  of  thermocouple  outputs  that  could  be  recorded  was 
limited  to  four  by  the  available  number  of  tape  recorder  channels.   A 
maximum  of  seven  channels  can  be  installed  on  the  recorder  used  (TMC 
Magnetic  Tape  Recording  System,  Model  700,  NE  1346)  but  only  four  plug-in 
circuits  were  available.   The  tape  recorder  can  be  operated  at  four 
different  tape  speeds,  1  7/8,  3  3/4,  15  and  30  inches  per  second. 

4.3.2  Analyzing  Equipment 

The  tape  recorder  and  the  analog  computer  were  used  in  the  record 
analysis  in  conjunction  with  a  hybrid  correlation  computer  (TMC  Correlation 
Computer,  Model  257,  NE  1363  with  TMC  Computer  of  Average  Transients, 
CAT  400C,  NE  1364).  A  photograph  of  the  arrangement  is  presented  in  Fig. 
21. 

The  thermocouple  outputs  were  analyzed  in  sets  of  two  on  the 
correlation  computer  at  a  higher  tape  speed  than  that  used  during  recording. 
Briefly,  the  functioning  of  the  hybrid,  analog-digital,  correlation  compu- 
ter is  as  follows.   One  of  the  inputs,  which  have  to  be  in  the  range  +  3  V, 
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is  digitized  and  the  value  stored  in  a  256  location  ferritic  core  memory, 
where  up  to  255  values  of  the  same  input  at  past  times  spaced  by  the 
delay  increment  time  are  kept.   During  each  delay  increment  time  the 
whole  memory  of  the  past  values  of  the  first  input  is  scanned  and  each 
value  is  multiplied  by  the  instantaneous  value  of  the  other  input.   The 
products  thus  obtained  are  accumulated  in  the  memory  of  the  Correlator 
of  Average  Transients  (CAT) .   In  the  first  location  of  the  CAT  memory  the 
products  of  the  simultaneous  values  of  the  two  inputs  are  accumulated, 
in  the  second  location  the  products  of  the  values  spaced  by  one  delay 
increment,  and  so  on.   The  precision  of  the  digitalization  of  the  inputs 
is  claimed  by  the  manufacturer  to  be  one  part  in  512  of  the  3  V  maximal 
input.  A  discussion  of  the  correlator  response  is  given  in  Appendix  C, 
and  an  outline  of  the  correlation  technique  in  Appendix  B. 

Three  output  devices  are  available  for  the  correlation  computer: 
punched  paper  tape,  analog  output  and  teletype  output.   A  Tally  Co. 
paper  tape  punch  was  normally  used  and  a  plot  of  the  correlation 
function  was  produced  on  an  X-Y  plotter,  mainly  for  identification  purposes. 
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5.   EXPERIMENTAL  PROCEDURE 

The  experimental  procedure  for  the  two  measurements,  the  thermal 
diffusivity  measurement  and  the  contact  resistance  measurement,  was 
identical  except  for  the  specimen  manipulations. 

A  special  device  was  used  for  the  alignment  of  the  thermal 
contact  resistance  specimens.   Two  polyethylene  half  rings,  3/4  inch 
wide  were  placed  around  the  contacting  sections  of  the  two  specimens. 
The  ring  was  tightened  hy  a  steel  band  with  a  sliding  latch.   When 
the  locking  piece  was  in  place  the  two  specimens  were  tightly  connected 
and  could  be  manipulated  as  a  unit. 

To  install  the  specimens  the  loading  lever,  the  top  plate,  and  the 
cylindrical  part  of  the  experiment  vessel  were  removed,  as  well  as  the 
top  cooler.   The  specimen  rod  was  pressed  into  the  bottom  cooler  and 
the  top  heater  and  cooler  were  installed.   The  specimen  thermocouple 
wires  were  soldered  to  permanent  extension  leads  in  the  bottom  compart- 
ment of  the  vessel.   The  cylinder  and  the  top  plate  were  assembled  with 
the  loading  ram  resting  on  the  top  cooler.   When  the  loading  lever  was 
replaced  and  the  specimens  were  secured  in  position  by  force,  the  sliding 
latch  was  released  by  pulling  a  wire.   The  parts  of  the  aligning  ring 
thereafter  were  not  in  contact  with  the  specimens.   The  desired  load  was 
applied  and  the  heaters,  coolers  and  the  microswitch  cam  drive  were  set 
into  operation. 

The  temperatures  of  both  the  guard  sleeve  and  the  specimen(s)  were 
measured  with  the  potentiometer.   The  guard  sleeve  heater  and  cooler 
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were  adjusted  so  that  the  guard  sleeve  temperature  matched  the  mean 
specimen  temperature  at  the  same  vertical  position  as  closely  as 
possible.   For  specimens  of  different  conductivity,  the  slope  could 
not  be  matched  completely,  but  the  differences  were  kept  small  (Fig.  22). 
When  there  was  no  more  temperature  drift,  the  specimen  thermocouples 
were  connected  to  the  amplifying  and  filtering  circuits.   The  bias 

voltages  were  adjusted  and  appropriate  amplification  chosen.   The 

3 

amplification  varied  from  10  for  the  thermocouple  closest  to  the 

oscillation  source  to  10  for  the  distant  thermocouples.  When  very  high 
amplifications  were  used  it  was  difficult  to  retain  the  amplified  signal 
within  the  1  volt  input  range  required  by  the  tape  recorder,  as  a  devia- 
tion of  ,1°C  would  drive  the  signal  beyond  that  limit.  Most  of  the  time 

it 

amplifications  of  up  to  2x10  were  used. 

Next  the  recorder  was  started,  using  a  recording  speed  of  1  7/8 
inches  per  second. 

Records  of  2  to  2  1/2  hours  duration  were  usually  taken,  which 
represents  80  to  100  periods  of  the  2.5  minute  period  waves  (the  longest 
period  used) .   The  recording  was  mostly  done  at  night  or  on  weekends  when 
the  disturbances  in  the  water  and  power  supply  were  small.   During  busy 
switching  hours  the  amplified  thermocouple  output  would  regularly  drift 
out  of  range  due  to  the  change  of  line  voltage  and/or  change  of  water 
temperature  and  pressure.   In  Fig.  23  typical  temperature  records 
obtained  are  shown. 

The  thermocouple  number  1  record  contains  higher  harmonic  components. 
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5°C 


Thermocouple   #1 


1°C 


Thermocouple   #2 


.1°C 


Thermocouple  #3 


1  min 


Figure  23.   Temperature  records. 
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The  ratio  of  the  next  higher  harmonic  to  the  fundamental  mode  is  estimated 
at  1/.172  from  the  autocorrelation  function.   No  significant  higher 
harmonic  was  found  in  the  thermocouple  number  2  record.   There  is  a 
significant  noise  component  in  the  output  of  the  thermocouple  number  4 
record.  Note  the  different  amplifications. 

The  magnetic  tape  was  later  played  back  at  16  times  the  recording 
speed  (i.e.  at  30  inches  per  second)  and  the  individual  records  or  pairs 
of  records  were  analyzed  on  the  correlation  computer.   The  autocorrelations, 
which  should  provide  the  amplitude  information,  were  later  omitted  as  it 
was  found  that  the  amplitude  response  of  the  correlator  is  quite  unre- 
liable. The  correlation  function  was  subsequently  punched  on  the  paper 
tape  and  an  analog  output  was  obtained  (an  example  is  presented  in 
Fig.  34) .  The  paper  tape  was  converted  to  magnetic  tape  and  subsequently 
to  punched  cards.   These  cards  in  turn  served  as  the  input  for  the 
correlation  function  analysis  code  described  in  Appendix  F. 
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6.  DATA  ANALYSIS 

6.1  Analysis  of  Temperature  Records 

The  purpose  of  the  data  collecting  system  in  the  temperature  wave 
experiment  is  to  obtain  the  information  on  the  temperatures  at  various 
positions  along  the  specimen  either  continuously  in  time  or  at  certain 
time  intervals.   In  this  investigation  the  temperatures,  represented  by 
the  thermocouple  voltage  outputs,  were  recorded  continuously  on  a  magnetic 
tape  recorder.   From  such  simultaneous  records  of  thermocouple  outputs  the 
information  on  the  phase  difference  and  the  amplitude  ratio  between  them 
is  to  be  obtained.  Various  methods  can  be  employed  for  this  purpose.  At 
present  a  relatively  convenient  method  of  analyzing  voltage  records  con- 
taining noise  components  (and  therefore  requiring  averaging  over  long 
periods  of  time)  is  available  which  utilizes  correlation  computers  with 
analog  input.   Previously,  the  correlation  technique,  the  theory  of  which 
is  briefly  outlined  in  Appendix  B,  was  possible  only  for  digitized  data. 

The  cross  correlation  function  is  defined  as  the  average  product  of 
two  records  with  the  delay  time  between  them  as  the  argument 


-    t  . 
$   (t)  -  lim  ~|  fjt   f(t)g(t+x)dt  (37) 
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where  f(t)  and  g(t)  are  the  two  records  and  t  the  delay  time  for  the 
record  g(t).  A  record  can  also  be  correlated  with  itself,  in  which  case 
the  above  function  is  termed  the  autocorrelation  function. 

In  Appendix  B  the  important  relationship  between  the  Fourier 
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transform  of  the  correlation  function  and  the  original  record  transforms 
are  derived, 


H*ff(T)}  =  ^f^{f(t»|?-  (38) 

i     i 


>{«f8(T>}  -  li^~F{i(t)}-'F{g(t)}  (39) 


where  the  symbol  F  denotes  the  Fourier  transform  and  the  asterisk  the 
complex  conjugate  value.   Thus  from  the  correlation  functions  (auto- 
and  cross-)  the  amplitude  and  phase  relationship  of  the  frequency 
components  of  the  records  are  available. 

The  correlation  of  the  records  with  one  prominent  frequency,  as 
are  obtained  in  the  temperature  wave  measurements,  is  a  special  case 
of  the  general  correlation  and  the  pertaining  relations  can  be  derived 
directly,  as  is  shown  in  Appendix  B,  Eq.  (B  13)  to  Eq.  (B  17).   The  lag 
between  the  harmonic  functions  is  in  this  case  retained  as  the  lag  of 
the  cosine. 

For  finite  and  discrete  correlation  functions  certain  limitations 
have  to  be  considered.   The  Fourier  transform  is  not  applicable  but  is 
reduced  to  the  Fourier  series  representation.   The  fact  that  the  fre- 
quencies in  the  Fourier  series  are  spaced  by  1/T,  where  T  is  the  length 
of  the  correlation  function,  reflects  the  uncertainty  with  which  the 
frequency  can  be  determined.   The  discrete  nature  of  the  correlation 
function  makes  it  impossible  to  distinguish  between  the  basic  frequency 
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and  its  aliases,  2kf  -f,  2kf  +f ,  where  k  is  any  integer  and  f  is  the 
Nyquist  frequency,  fN  =  1/2At .   (Ax  is  the  distance  between  the  points). 
This  effect  can  be  predicted  from  the  development  of  the  Fourier  series 
on  a  set  of  discrete  points,  but  the  development  is  not  presented  here 
since,  for  the  correlations  performed,  the  Nyquist  frequency  was  in  the 
range  where  the  filtering  circuits  eliminated  any  significant  component 
from  the  records. 

The  information  that  the  record  and  consequently  the  correlation 
function  contains  a  discrete  prominent  frequency,  led  to  a  different 
frequency  analysis  approach  than  the  Fourier  series  expansion.   A  linear 
combination  of  a  polynomial  and  one  or  more  harmonic  terms  were  fitted 
to  the  correlation  function  in  the  least  sum  of  squares  sense.   The 
polynomial  part  was  introduced  to  represent  the  decaying  exponential 
component  in  the  case  of  the  autocorrelation  function  of  a  noisy  record. 
This  model  was  found  to  be  correct  as  the  residuals  were  random  in  the 
great  majority  of  cases.   The  frequency  so  determined  should  still  be 
considered  uncertain  in  the  interval  Af  =  1/T  where  T  is  the  length  of 
the  correlation  function,  in  the  sense  that  should  another  harmonic  be 
present  within  the  interval  it  would  not  be  possible  to  determine  them 
simultaneously.   In  this  case  the  one  harmonic  model  would  prove  incorrect 
(many  harmonics  would  have  to  be  introduced) ,  and  the  coefficients  of  the 
two  harmonics  would  be  highly  correlated. 

An  anomaly  was  discovered  in  the  response  of  individual  channels 
of  the  correlation  computer  and  an  additional  function  was  introduced  as 
a  component  of  the  linear  combination  to  be  fitted  to  the  correlation 
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function.   The  response  of  the  correlator  is  discussed  in  Appendix  C 
and  a  plot  of  the  correlator  background  function  is  presented  in  Fig.  24. 
The  inclusion  of  this  background  function  reduced  the  number  of  cases 
in  which  the  residuals  were  not  random  and  a  further  refinement  might 
completely  eliminate  the  occurrence  of  such  cases.   In  this  regard,  it 
is  expected  that  the  individual  channel  response  is  also  nonlinear  as  is 
the  case  with  the  overall  amplitude  response  (cf.  Appendix  C) .  Although 
the  peaks  of  the  background  function  were  always  found  at  the  same  location 
their  relative  amplitude  may  not  be  the  same  for  different  inputs. 

The  least  squares  fit  of  the  analytic  form  to  the  correlation  function 
provided  the  amplitude  and  phase  of  the  prominent  harmonic  in  the  correlation 
function.   The  thermal  diffusivity  and  the  thermal  contact  resistance  were 
calculated  from  the  phase  angle  information,  as  described  in  the  following 
section.   The  amplitude  information  was  not  used  as  it  was  considered  too 
unreliable  in  view  of  the  findings  on  the  correlator  response  (Appendix  C) . 

The  theory  of  the  regression  analysis  (least  square  fitting)  is  given 
in  Appendix  E  and  the  computer  program  that  performs  the  task  is  described 
in  Appendix  F. 

6.2  Thermal  Diffusivity  Determination 

The  thermal  diffusivity  specimens  were  of  sufficient  length, 
compared  with  the  attenuation  distance  for  the  periods  used,  that  the 
no-reflection  approximation  for  the  behavior  of  the  temperature  wave 
could  be  used  (Eq.    (13)).   For  the  thermal  diffusivity   calculation 
the  following  expression  is  used 
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_  a    ,A1*2    n  ,41  2  .... 

a  "  2  W     "to  (A*>  (40) 


where  Al  is  the  distance  between  the  thermocouples  and  A d>  is  the  phase 
lag  between  the  thermocouple  outputs. 

6.3  Contact  Conductance  Determination 

The  results  derived  in  section  2.2  were  used  for  the  determination 
of  the  contact  conductance  parameter  from  the  phase  lag  information.   The 
assumption  of  a  discontinuity  in  an  infinite  medium  is  valid  as  the 
distance  from  the  boundaries  was  several  times  greater  than  the  attenuation 
length. 

The  contact  conductance  h  cannot  be  determined  solely  from  the  phase 
or  amplitude  ratio  and  the  distances  between  the  measuring  points,  which 
is  evident  also  from  its  dimensionality.   Only  the  dimensionless  parameter 
Y  =  h/kX  can  be  obtained,  and  the  knowledge  of  the  thermal  conductivity 
and  thermal  diffusivity  (contained  in  A  together  with  the  period  of 
oscillations)  is  needed.   Furthermore  for  the  determination  of  y  the 
knowledge  of  the  material  properties  ratio  r  (Eq.  (36))  is  necessary. 
This  latter  quantity  can  be  obtained  from  the  ratio  of  the  thermal  conduct- 
ivities and  the  ratio  of  thermal  diffusivities . 

From  the  results  of  the  present  experiments  the  value  of  y   was 
obtained  and  the  thermal  contact  conductance  was  calculated  using  the 
values  for  conductivity  from  the  literature.   The  ratio  of  the  conduct- 
ivities of  the  two  materials  was  estimated  from  the  temperature  gradient 
ratio  in  the  two  specimens  (since  the  same  amount  of  heat  flows  through 
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both) . 

The  thermal  contact  conductance  parameter  was  determined  from  the 
phase  difference  between  the  two  thermocouples  closest  to  the  discon- 
tinuity, which  is  most  sensitive  to  this  parameter.  As  no  explicit 
expression  is  available  the  parameter  was  determined  by  iteration  in 
a  computer  code.  The  input  data  were  the  thermal  diffusivities  of  the 
two  materials,  the  ratio  of  the  conductivities,  the  distances  of  the 
thermocouples  from  the  discontinuity  and  the  phase  difference.   For 
each  value  of  y   the  reflected  and  transmitted  wave  amplitude  and  phase 
at  the  discontinuity  were  calculated  (system  of  equations  Eqs.  (32)) 
and  from  these  the  phase  difference  at  the  two  thermocouple  positions. 
The  procedure  was  iterated  until  a  satisfactory  convergence  was  reached 
(the  difference  of  the  two  consecutive  estimates  less  than  1%). 
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7.   EXPERIMENTAL  RESULTS 

7.1  Thermal  Diffusivity  Measurements 

The  thermal  diffusivity  of  an  aluminum  alloy,  2024-T4,  and  Armco 
iron  was  measured  using  the  temperature  wave  method  described  herein. 
The  data  on  these  materials  and  the  experimental  results  are  given 
in  Table  II.   Values  for  the  thermal  diffusivity  computed  from  data 
obtained  from  the  literature  are  included  for  comparison.   The  experi- 
mental values  agree  very  well  with  the  literature  data. 

7.2  Thermal  Contact  Conductance  Measurements 

The  thermal  contact  conductance  of  aluminum-aluminum  and  aluminum- 
iron  contacts  was  measured.   The  measurements  were  performed  at  different 
applied  pressures.   Four  specimen  pairs  were  used  and  12  complete  data 
recordings  and  evaluations  were  performed. 

In  the  calculations  of  the  contact  conductar.ee  parameter  the  values 

of  the  thermal  diffusivity  measured  as  a  part  of  this  investigation  were 

2  2 

used,  the  values  being  .485  cm  /sec  for  aluminum  and  .208  cm  /sec  for 

iron.   The  ratio  of  the  thermal  conductivities  of  the  two  materials  was 

estimated  from  the  temperature  gradients  in  the  two  specimens  (cf.  Fig.  22), 

but  as  the  experiment  was  not  designed  for  this  purpose  the  data  from  the 

literature  were  used.   From  the  literature  (1,  5)  the  ratio  is  1.65  while 

the  mean  of  the  graphical  estimates  was  1.69. 

The  information  on  the  specimen  pairs  used  and  the  results  of  the 

correlation  runs  are  summarized  in  Tables  III  to  VI.   The  thermocouple 

position  1  is  closest  to  the  oscillating  heater,  the  thermocouple  3  last 
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Table  II 
Thermal  diffusivity  measurement  results 

Specimen  material:  Aluminum  2024-14.   Oscillation  period:  1.92  min. 

Thermocouple        Thermocouple        Phase  difference       1    , 

positions  distance  ,      ..    .  ,  . 

r  ,  .  (radians;  (cm) 

(cm) 

1-2  6.33  1.5107  4.20 

1-3  12.8  3.0095  4.26 

1-4  19.0  4.5528  4.17 

Average  L1  .  =  4.21 


2 

Thermal  diffusivity  measured   .485  cm  /sec        . 

Thermal  diffusivity  from  literature  data*   .465  cm  /sec 


Specimen  material:  Armco  Iron.   Oscillation  period:   3.25  min. 

Thermocouple        Thermocouple        Phase  difference  L1  , 

positions           distance             ,   ,.    .  ,  , 

*                                                 ,     ,               (radians)  (cm) 
(cm) 

1-2              6.40              1.809  3.55 

1-3             12.80              3.454  3.71 

1-4             19.2               5.403  3.56 

Average**  L  ,  =3.55 


1/e 


Thermal  diffusivity  measured   ,203  cm  /sec        „ 
Thermal  diffusivity  from  literature  data*   .208  cm  /sec 


3 

*Aluminum  (1,  5):   density  =  2.73  g/cm  ,  specific  heat  =  .23  cal/g°C 

thermal  conductivity  =  .29  cal/cm  sec°C 

3 
Armco  Iron:   density  =  7.86  g/cm  ,  specific  heat  at  30°C   .108  cal/g 
thermal  conductivity  at  30°C  =  .176  cal/cm  sec  °C. 

^Thermocouple  at  position  number  3  was  found  defective  so  the  data 
for  this  point  were  excluded  from  evaluation. 
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Table  III 


Thermal  contact  conductance  measurements,  pair  A11-A12 


Specimen  1 

Specimen  2 

Al  1 

2024-T4 
77 

Al  2 

2024-T4 

77 

Identification  code 

Material 

Hardness,  HR^^  mQ   fcp) 

Test  surface  roughness,  pin  rms        15  10 

Thermocouple  no.  3  (in  specimen  1),  distance  from  interface  1.27  cm 
Thermocouple  no.  4  (in  specimen  2),  distance  from  interface  1.27  cm 


Load  250  lb  (pressure  320  psi),  parallel  lay 

Date  of  recording  12/7/1968,  Period  of  oscillation  1.92  min. 


Correlation,  thermocouple  no. 


Phase  difference,  radians 


1-2 

1.2057 

1-2 

1.2018 

1-3 

2.3294 

1-4 

3.0780 

2-4 

1.8510 

3-4 

.6410 

Average  phase  diffi 

2rence 

1  - 

-  3 

2.3318 

1  - 

3  - 

-  4 

.6718 

4  3.0345 


Load  580  lb  (pressure  790  psi) ,  parallel  lay 

Date  of  recording  12/8/1968,  Period  of  oscillation  1.92  min. 


Correlation,  thermocouple  no. 

1-2 
1-3 
1-3 
1-4 
2-4 
3-4 

Average  phase  difference 


Phase  difference,  radians 


1.1990 

2.3326 

2.3524 

2.9354 

1.8503 

0.7150 

1  - 

-  3 

2.3425 

1  - 

3  - 

-  4 

0.6531 

2.9955 
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Table  IV 
Thermal  contact  conductance  measurements,  pair  A3  4-A15 


Specimen  1 

Specimen  2 

Al  4 

Al  5 

2024-T4 

2024-T4 

78 

77 

Identification  code 

Material 

Hardness,  HR^^  m   fcp) 

Test  surface  roughness, Pin  rms         30  25 

Thermocouple  no.  3  (in  specimen  1),  distance  from  interface  1.27  cm 
Thermocouple  no.  4  (in  specimen  2),  distance  from  interface  1.32  cm 

Load  300  lb  (pressure  380  psi) ,  parallel  lay 

Date  of  recording  12/9/1968,  Period  of  oscillation  1.92  min. 

Correlation,  thermocouple  no.     Phase  difference,  radians 

1-2 
1-3 
1-3 
1-4 
1-4 
1-4 
3-4 

Average  phase  difference  1-3   2.2734     1-4   3.0655 


Load  630  lb  (pressure  850  psi) ,  parallel  lay 

Date  of  recording  12/8/1968,  Period  of  oscillation  1.92  min. 

Correlation,  thermocouple  no.     Phase  difference,  radians 

1-2 
1-3 
1-3 

1-4 
1-4 
3-4 

Average  phase  difference   1-3   2.2124    1-4   2.9973 


1.1743 

2.2673 

2.2579 

3.1099 

2.9811 

3.0640 

0.8344 

1  - 

3 

2.2734 

1  - 

3  - 

4 

.7913 

1.1796 

2.2518 

2.1731 

2.9892 

3.0143 

0.7712 

1  - 

•  3 

2.2124    1  - 

3  • 

■  4 

.7849 
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Table  IV.  continued 

Load  1750  lb  (pressure  2250  psi),  parallel  lay 

Date  of  recording  12/10/1968,  Period  of  oscillation  1.92  min 

Correlation,  thermocouple  no.  Phase  difference,  radians 

1  -  2  1.3596 

1-3  2.2045 

1  -  A  2.8784 

1-4  2.9026 

3-4  0.6703 

Average  phase  difference   1-3  2,2045     1-4   2.8853 

3-4  .6808 
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Table  V 

Thermal  contact  conductance  measurements,  pair  A13-I4 

Specimen  1     Specimen  2 

Identification  code  Al  3  14 

Material  2024-T4       Armco-Iron 

Hardness,  fflj^, _  1QQ   fcp)  76  76 

Test  surface  roughness , uin  rms         50  35 

Thermocouple  no.  3  (in  specimen  1),  distance  from  interface  1.27  cm 
Thermocouple  no.  4  (in  specimen  2),  distance  from  interface  1.29  cm 

Load  250  lb  (pressure  320  psi) ,  cross  lay 

Date  of  recording  1/2/1969,  Period  of  oscillation  2.5  min. 

Correlation,  thermocouple  no.     Phase  difference,  radians 

1-2  1.0685 

1-3  2.1000 

1-4  2.8505 

1-4  2.7847 

2-3  1.0664 

2-4  1.8270 

3-4  0.8180 

Average  phase  difference   1-3   2.1174    1-4   2.8875 

3-4   0.7701 


Load  550  lb  (pressure  750  psi) , 

Date  of  recording  1/2/1969,  Period  of  oscillation  2.5  min. 

Correlation,  thermocouple  no.  Phase  difference,  radians 

1-2  1.0529 

1-3  2.0550 

1-4  2.7492 

2-3  1.0163 

2-4  1.7032 

3-4  0.6782 

Average  phase  difference   1-3  2.0616    1-4   2.7480 

3-4  .6864 
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Table  V.  continued 


Load  900  lb  (pressure  1170  psi) ,  cross  lay 

Date  of  recording  1/2/1969,  Period  of  oscillation  2.5  min 

Correlation,  thermocouple  no.  Phase  difference,  radians 

1-2  1.0719 

1-3  2.0912 

1  -  4  2.7787 

2-3  1.0413 

2-4  1,7324 

3-4  0.6757 

Average  phase  difference   1-3  2.1022    1-4   2.7867 

3-4  .6844 
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Table  VI 
Thermal  contact  conductance  measurements,  pair  A14-I1 


Specimen  1 

Specimen  2 

Al  4 

I  1 

2024-T4 

Armco  Iron 

78 

72 

Identification  code 

Material 

Hardness,  H^^,,  m   ^ 

Test  surface  roughness , uin  rms  15  10 

Thermocouple  no.  3  (in  specimen  1),  distance  from  interface  1.27  cm 
Thermocouple  no.  4  (in  specimen  2),  distance  from  interface  1.315  cm 

Load  250  lb  (pressure  320  psi) ,  cross  lay 

Date  of  recording  1/4/1969,  Period  of  oscillation  2.5  min. 

Correlation,  thermocouple  no.     Phase  difference,  radians 

1-2  1.0573 

1-3  2.0649 

1-4  2.7313 

2-3  1.0097 

2-4  1.6963 

3-4  0.6773 

Average  phase  difference   1-3   2.0665    1-4   2.7425 

3-4   0.6761 


Load  550  lb  (pressure  700  psi),  cross  lay 

Date  of  recording  1/5/1969,  Period  of  oscillation  2.5  min. 

Correlation,  thermocouple  no.     Phase  difference,  radians 

1-2  1.0912 

1-3  2.0706 

1-4  2.7343 

2-3  1.0124 

2-4  1.6706 

3-4  0.6736 

Average  phase  difference   1-3   2.0866    1-4   2.7520 

3-4   0.6654 
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"Table  VI.  continued 


Load  900  lb  (pressure  1140  psi) ,  cross  lay 

Date  of  recording  1/4/1969,  Period  of  oscillation  2.5  min. 


Correlation,  thermocouple  no. 

1-2 
1-3 
1-4 
2-3 
2-4 
3-4 


Phase  difference,  radians 

1.0611 
2.0599 
2.7214 
1.0134 
1.6594 
0.6629 


Average  phase  difference 


2.0671 
0.659S 


2.7270 


Load  1200  lb  (pressure  1510  psi), 

Date  of  recording  1/4/1969,  Period  of  oscillation  2.5  min. 


Correlation,  thermocouple  no. 

1-2 
1-3 
1-4 
2-3 
3-4 

Average  phase  difference   1 

3 


Phase  difference,  radians 

1.0290 
2.0132 
2.7300 
0.9735 
0.6974 


2.0079 
0.7098 


2.7176 
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in  specimen  1  before  the  interface  and  thermocouple  4  is  in  specimen  2. 

The  resulting  thermal  contact  conductance  values  are  presented  in 
Table  VII  and  in  Fig.  25.   In  Fig.  25  the  experimental  values  and 
theoretical  prediction  as  given  in  reference  (15)  are  plotted  for  the 
aluminum  2024-T4  -  Armco  iron  contact.   The  data  and  the  theoretical 
curve  are  for  a  rms  roughness  of  150 u  in,  which  is  much  higher  than 
that  of  the  specimens  used  in  this  work. 

The  phase  difference  between  the  two  thermocouples  closest  to  the 
interface,  used  in  the  calculation  of  the  thermal  conductance  parameter, 
is  the  average  value  for  this  phase  difference  as  indicated  by  the 
difference  of  phase  between  correlations  1-3  and  1-4;  2-3  and  2-4  and  the 
direct  correlation  3-4.   The  correlations  with  the  thermocouples  1  and  2 
give  a  more  reliable  estimate  than  the  direct  correlation  of  the  two 
thermocouples  close  to  the  interface  (3  -  4)  alone,  as  the  records  1  and 
2  have  much  smaller  relative  noise  and  provide  good  reference  records. 
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Table  VII 


Summary  of  thermal  contact  conductance  results 


Specimen 

Applied 

Conductance 

Thermal 

contact 

pair 

pressure 

parameter 

conductance 

(psi) 

(Btu/hr  ft2 
X10 

°F) 

(cal/cm  se 

Al  1-A1  2 

320 

6.8 

6.4 

.47 

Tl        1! 

740 

9.5 

8.8 

.65 

Al  4-A1  5 

380 

2.3 

2.16 

.159 

810 

2.4 

2.24 

.165 

2250 

7.6 

7.1 

.525 

Al  4-1  1 

320 

10.0 

8.15 

.60 

710 

12.8 

10.5 

.77 

1150 

14.9 

12.1 

.89 

1540 

5.9 

4.7 

.355 

Al  3-1  4 

320 

3.2 

2.8 

.206 

710 

10.7 

8.7 

.64 

1150 

11.3 

9.2 

.68 

°c) 
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Figure. 25.   Contact  thermal  conductance  experiment  results. 
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8.   DISCUSSION  OF  RESULTS  AND  EXPERIMENT 

The  experiments  performed  demonstrated  the  feasibility  of  deter- 
mining the  thermal  contact  conductance  by  the  temperature  wave  method, 
the  determination  of  the  thermal  diffusivity  having  been  accomplished 
before.   In  the  following  the  interpretation  of  results  is  given  and 
some  shortcomings  of  the  present  work  are  discussed. 

The  thermal  diffusivity  measurement  results  agree  very  well  with 
the  published  data  for  the  materials  examined.   The  choice  of  the 
aluminum  alloy  in  the  tempered  condition  was  not  the  most  appropriate 
for  this  purpose  as  Its  properties  are  given  only  as  approximate  while 
for  the  alloy  in  the  annealed  condition  more  exact  data  are  available. 
The  difference  in  results  is  within  5%  and  the  agreement  can  be  considered 
as  very  good.   The  measured  value  of  the  thermal  diffusivity  was  used  in 
the  further  analysis  as  it  is  considered  more  applicable  for  the  particular 
material  than  the  general  value  from  the  literature. 

A  minor  deficiency  of  the  experimental  procedure  may  be  illustrated 
with  the  data  obtained  in  the  thermal  diffusivity  measurement  of  aluminum. 
Due  to  the  necessary  connection  of  one  wire  of  each  thermocouple  to  the 
analog  computer  ground,  additional  loops  were  created  between  the  thermo- 
couples as  they  were  mostly  in  good  electrical  contact  with  the  specimen 
(loop  ABC  in  Fig.  26c).   In  these  loops  current  flowed  due  to  the  tempera- 
ture difference  between  the  two  thermocouple  positions,  where  the  thermo- 
couple copper-specimen  material  was  created.   The  minute  current  caused 
a  voltage  drop  in  each  leg  CAB  and  AC),  which  was  added  or  subtracted 
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Figure  26. a.   Deviations  of  phase  angle  from  linearity 


Figure  26. b.   Phase  angle  and 
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thermocouple  outputs 


AMPLIFYING  CIRC. 


i 


REFERENCE 
JUNCTION 


Figure  26. c.   Thermocouple  wire  connections. 
Figure  26.   Thermocouple  interdependence. 
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from  the  voltage  generated  by  the  respective  thermocouples.   This  was 
true  as  well  for  the  steady  part  as  for  the  oscillatory  part  of  the 
temperature.  Among  the  oscillatory  voltages  the  dominant  was  the- 
output  of  the  thermocouple  with  the  highest  temperature  amplitude  and 
it  influenced  the  others,  though  this  effect  was  small.   The  phase 
relation  of  the  thermocouples  in  the  aluminum  thermal  diffusivity 
specimen  as  well  as  the  effect  on  the  phase  deviation  from  the  mean 
phase  slope  can  be  seen  from  Fig.  26  a  and  b.   The  effect  is  in  this 
case  hardly  significant,  but  was  noticed  in  other  results  also.   The 
uncertainty  intervals  are  the  confidence  limits  for  the  phase  determin- 
ation from  the  correlation  function  (cf.  Appendix  E) .   Though  small 
this  effect  should  be  avoided  by  securing  insulation  of  the  thermocouples 
from  the  specimen.  The  magnitude  of  the  deviations  was  still  smaller 
for  the  thermal  contact  conductance  measurements  as  the  ratios  of 
amplitude  were  more  modest. 

The  thermal  contact  conductance  parameter  measurements  rely  on  the 
ability  to  determine  the  deviations  from  linearity  of  the  phase  relations 
and  are  thus  less  accurate  than  the  diffusivity  measurements.   In  Fig. 
25  the  measured  contact  conductances  are  plotted  together  with  the 
measurements  reported  in  reference  (15),  which  were  the  only  measurements 
found  in  the  literature  for  the  material  pair  used  and  with  air  as  the 
interstitial  medium.   It  should  be  noted  that  the  contacting  surfaces 
for  which  the  contact  conductances  are  reported  in  reference  (15)  were 
much  rougher  than  the  ones  used  here,  the  rms  deviation  being  150  pin. 
compared  to  10  to  50  pin.  used  in  this  investigation.   For  surfaces  of 


79 


different  roughness  the  difference  of  the  contact  conductance  may  amount 
to  up  to  one  order  of  magnitude,  as  reported  in  reference  (24). 

A  deficiency  of  the  thermal  contact  conductance  measurement  program 
was  that  no  detailed  information  on  the  contacting  surface  waviness 
was  available.   The  surfaces  were  machined  very  carefully  and  no  waviness 
could  be  detected  by  measurements  with  a  dial  indicator  at  several  points 
on  the  surface.  No  satisfactory  assurance  is  though  given  that  the 
surfaces  were  indeed  flat.   This  fact  does  not  imply  that  the  results 
for  the  thermal  contact  conductance  are  inaccurate  but  only  that  the 
comparison  with  other  experimental  results  is  hindered. 

A  few  experimental  points  evidently  fall  beyond  the  random  error 
of  the  measurements.   These  are  the  points  for  Al  3  -  I  4,  320  psi 
and  Al  4  -  I  1,  1540  psi.   The  documents  of  the  experiments  were  reviewed 
and  it  was  found  that  the  erroneous  points  were  obtained  from  records 
with  considerable  drift  while  in  the  other  records  the  drift  was  negli- 
gible.  It  should  be  noted  that  the  expressions  of  the  correlation 
theory  (Appendix  E)  are  derived  for  stationary  records  and  any  drift  of 
important  slope  will  cause  disturbance  in  the  phase  relationships  obtained. 
Care  should  be  therefore  taken  that  the  steady  condition  is  indeed 
obtained  and  maintained.   The  records  should  be  inspected  for  drift  and 
the  experiment  repeated  to  avoid  the  rather  uncomfortable  posterior 
rejection  of  erroneous  points  on  the  basis  that  they  do  not  agree  with 
other  results.   One  of  the  erroneous  points  (Al  3-14,  320  psi)  was 
obtained  in  an  experiment  where  the  matching  of  the  specimen  and  guard 


sleeve  temperatures  was  poor  (difference  of  2  to  3°C)  and  the  record 
also  contained  some  unstationarity ,  the  latter  probably  being  the  main 
source  of  error. 

The  experimental  points  show  the  expected  trend  of  higher  conduct- 
ance for  smoother  surfaces  and  higher  pressure.  Also  lower  slope  in 
the  conductance  -  pressure  diagram  can  be  expected  for  smoother  surfaces, 
as  the  measurements  indicate. 

The  thermal  contact  conductance  of  all-aluminum  pairs  also  falls 
in  the  range  of  aluminum-iron  results,  as  the  hardness,  roughness  and 
the  method  of  machining  the  surfaces  was  close.  The  choice  of  the 
parallel  lay  position  is  somewhat  inappropriate  as  in  this  case  the 
contacting  points  are  distributed  in  a  quite  different  pattern  than  in 
the  cross-lay  position. 

Evaluating  the  results  of  the  thermal  contact  conductance  measure- 
ments it  can  be  stated  that  satisfactory  accuracy  can  be  achieved  with 
this  method  if  sufficient  care  is  taken  to  maintain  stationary  conditions. 
The  ultimate  and  not  very  remote  limitation  to  accuracy  may  be  in  the 
present  experimental  equipment,  particularly  in  the  method  of  analyzing 
the  records  with  the  correlation  computer.   The  adequacy  of  the  method 
for  the  thermal  contact  conductance  measurements  is  further  supported 
by  the  fact  that  the  reproducibility  of  the  phenomenon  itself  is  on  the 
order  of  20  to  30%  error,  as  estimated  from  various  compilations  of  data 
(10,  15.). 


81 


9.   SCOPE  FOR  FURTHER  WORK 

The  method  of  thermal  contact  conductance  measurement  can  be  con- 
sidered established  and  the  investigation  of  the  phenomenon  itself 
should  be  pursued.   The  method  of  temperature  waves  may  be  particularly 
valuable  in  providing  new  insight  in  the  directional  effect  in  contact 
conductance.   In  this  respect  a  modification  of  the  experimental  apparatus 
should  be  considered,  which  would  allow  steady  state  and  temperature  wave 
measurements,  with  heat  flow  in  either  direction,  obtainable  without 
displacing  the  specimens. 

The  present  temperature  record  analysis  system  was  found  inadequate 
for  determining  the  amplitude  relations  of  the  records.   A  different 
kind  of  data  acquisition  system  might  be  considered,  or  the  present  one 
improved.   The  sampling  interval  of  the  temperatures  is  several  seconds 
(in  the  present  equipment  2.5  seconds)  and  the  records  are  not  excessive 
in  length,  so  that  the  digitizing  of  the  thermocouple  voltage  output  may 
be  appropriate  and  all  further  analysis  performed  in  the  digital  form. 
A  precision  analog-digital  converter  would  be  needed  (digital  voltmeter) 
and  a  digital,  computer  readable,  storage  (magnetic  tape). 

Such  arrangement  could  provide  for  more  measuring  points  in  the 
specimens  and  monitoring  of  the  temperatures  outside  the  specimen.   An 
additional  advantage  of  such  system  would  be  that  the  sampling  could  be 
synchronized  with  the  oscillation  generating  devices  (heater  and  cooler) 
and  so  an  ideal  reference  record  for  correlation  would  be  provided.   Such 
equipment  could  also  be  used  for  other  experiments  with  temperature  waves. 
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APPENDIX  A:  Transient  Response  of  Elemental  Heat  Channel 

In  the  investigation  of  the  phase  change  and  the  amplitude  of  the 
reflected  and  transmitted  waves  at  the  interface  of  two  solids  a  constant 
heat  transfer  coefficient  is  assumed  and  no  heat  capacity  is  assigned  to 
the  resistance.   The  resistance  is  due  to  the  constriction  of  the  heat 
flow  area  and  the  flow  lines  are  displaced  to  a  certain  distance  from 
the  actual  contact,   so  some  thermal  capacity  may  be  associated  with  it. 

In  order  to  obtain  the  step  input  response  of  an  elemental  channel 
with  constriction  resistance  the  shape  of  the  heat  flow  channel  is 
approximated  by  the  shape  shown  in  Fig.  27,  as  suggested  by  Heasley  (17). 
In  this  approximation  the  area  of  the  conical  part  increases  with  the 
square  of  the  distance  from  the  apex  of  the  cone.   It  is  not  implied 
that  the  geometrical  shape  of  the  actual  elemental  channel  is  similar 
to  the  one  presented  but  only  that  the  heat  is  assumed  to  flow 
through  a  path  that  can  be  represented  in  this  way.   The  transient 
response  of  this  elemental  channel  shall  be  compared  to  that  of  a 
uniform  channel  (rod)  with  the  step  rise  of  temperature  applied  to  one 
end  through  a  constant  heat  transfer  coefficient  (Fig.  27). 

For  the  constant  cross  section  channel  the  problem  is  formulated 
as  follows. 
The  governing  equation  is 


a2l<*,t)     1  w(,.t) 

3x2   "a    at  CA-1) 
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Figure  27. a.   Constricted  heat  channel  model 


Figure  27. b.   Constant  cross  section  heat  channel  model  with 
input  resistance. 


Figure  27.   Two  models  of  elemental  heat  channel. 
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where  a  =  — —  is  the  thermal  diffusivity. 
CPP 


The  boundary  condition  applied  at  x  =  0  is 


-k  3TfX't)  =  h(T  -  T(x,t))     x  =  0,  t  >  0  (A-2) 

3x        o 


where  T  is.  the  temperature  level  applied  at  t  =  0,  and  h  is  the 
heat  transfer  coefficient  at  the  surface. 
The  initial  condition  is 


T(x,0)  =  0 


and  the  other  boundary  condition 


lim  T(x,t)  =  0 


Applying  the  Laplace  transform  to  the  governing  equation  we  have 


^%^i  =  ^T(x,s) 
dx2 


(A-3) 


which  has  the  general  solution 


T(x,s)  =  C  e'a   +  C2 


/s         /s 
1=  X  +  Ce""'^ 


(A-A) 
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Since  the  temperature  should  be  bounded  for  large  x,  C.  =  0,  and 
application  of  the  boundary  condition  at  x  -   0  yields 


kf£c  e~'aX  =  h(T  /s  -  C,e"'a  X)  (A-5) 

I  a  2  o      2 


The 

transformed 

solution  is 

thus  found 

to 

be 

T(x 

s) 

hT 
o 

.*£  + 

la 
e 

h) 

(A-6) 


The  inverse  transform  of  this  expression  can  be  readily  obtained, 
however  this  is  not  true  for  the  next  case  and  thus  only  a  comparison 
of  the  transformed  solutions  shall  be  made. 

For  the  approximate  constricted  elemental  channel  (Fig.  27)  it 
shall  be  assumed  that  the  heat  flow  lines  in  the  conical  part  are 
straight  and  diverge  from  the  common  origin  on  the  axis  of  the  heat 
channel  at  x  =  0.   With  this  assumption  we  have  a  one  dimensional 
transient  problem  in  the  conical  region  also.   The  two  governing 
equations  are 


32V(*,Q  +  2  3V(x,t)  =  1  3V(x,t)     d  <  x  <  b  (A_7) 

.  .  2      x   3x      a   3t 
3x 


i>2T(*.t)  .  I  Wfct)  x  „  b  (A_8) 

3x2      °    3t 
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where  V  is  the  temperature  in  the  conical  region. 
The  boundary  and  the  initial  conditions  are 

T(x,0)  »  V(x,0)  =0  x  >  d 

V(d,t)  -  T  t  >  0 
o 

V(b,t)  =  T(b,t)  t  >  0 

|I(b>t)-f(b,t)  t>0 

lim  T(x,t)  =0  t  >  0                  (A-9) 
x-*» 


The  Laplace  transformed  equations  are 


1  d2(xV(x,s))  _  s  y(x    m   0    d<x<b  (A.10) 

x     ,2       a 
dx 


2- 
d  I(*'s)  -  -  T(x,s)  =0        b  <  x  (A-ll) 

dx       a 


The  solution  of  Eq.  (A-7)  satisfying  the  conditions  (A-9)  is 


T  d  r  ci         r 

V(x,s)  =  -5-  cosh([-  (x-d))  +  —  sinh([-  (x-d))         (A-12) 
sx       la         x       |a 


and  the  solution  of  Eq.  (A-8)  satisfying  the  conditions  Eq.  (A-9)  is 

T(x,s)  =  C2e  ,a  (A-13) 
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Using  the  boundary  conditions  at  x  ■  b  to  evaluate  the  constants  C. 
and  C„  we  obtain  for  the  temperature  in  the  uniform  part  of  the  rod 
(temperature  in  the  conical  part  being  of  lesser  interest) 

-f-  (x-b) 

l_de 
T(x,s)   = 


sb[sinh(P  (b-d))  +  cosh(]^  (b-d))   -  -±=  sinh(fi  (b-d)) 
I  cc  ia  |s^  i  a 


(A-14) 


To  compare  this  result  with  the  one  obtained  previously  we  should 
calculate  the  resistance  of  the  constriction  to  a  steady  heat  flow. 
Under  the  same  assumptions  as  stated  before  for  the  heat  flow  lines 
the  resistance  due  to  constriction  is 


R  .  b£p»  .  1/h  (a_i5) 


If  we  expand  the  hyperbolic  functions  in  Eq.  (A-14)  into  series  and 
retain  only  the  linear  terms  we  obtain 

T  (x,s)  =  <*><*-»  .    (A-16) 

s(I?  +  blb^dT) 

This  expression  is  identical  to  Eq.  (A-6)  if  the  stationary  constriction 
resistance  is  substituted  for  the  heat  transfer  resistance  1/h. 

In  order  that  the  truncation  of  the  series  be  acceptable  the  ratio 
of  the  lowest  order  term  neglected  to  the  terms  retained  should  be 
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much  less  than  one.   We  have  the  condition 


£<► 


d)  «  2 


s  <  ~  (A-17) 

(b-d)Z 


Thus  the  time  constant  of  the  transient  effect  is  of  the  order 


tt      .  ^  (A-18) 

trans.     2 


A  rough  estimate  of  typical  values  for  metallic  contacts  follows. 

2  2 

Using  the  ratio  of  apparent  to  real  area  of  contact  b  /d  =  100, 

2 

k  =  .2  cal/cm  sec  °C  and  h  =  .2  cal/cm  sec  °C,  the  expression  for 

the  steady  state  resistance  of  the  channel  model  with  contraction 
Eq.  (A- 15)  yields 


b  -  d  =  5— r  =  .1  cm 
h  b 


2 

Assuming  for  a  a  typical  value  for  metals,  .5  cm  /sec,  we  have 


t^      =  Order  (.1  sec) 
trans . 


Since  in  the  present  investigations  periods  of  oscillation  have  been 
employed  that  are  3  orders  of  magnitude  longer  than  this  time  it  is 
concluded  that  the  measured  contact  conductance  is  identical  to  the 
stationary  conductance  as  far  as  the  transient  effects  are  concerned. 
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APPENDIX  B:   Correlation  Technique. 

In  the  analysis  of  records  containing  random  noise  components,  the 
cross-  and  autocorrelation  techniques  have  proven  useful.   Kith  the 
development  of  convenient  hybrid  correlation  computers,  the  correlation 
analysis  of  a  much  wider  class  of  records  has  become  possible,  including 
high  frequencies  of  the  communication  technology,  whereas  before  only 
digitized  records  could  be  analyzed. 

B.l  Correlation  Function  Definition  and  Properties 

Although  actually  discrete  analysis  is  performed  on  records  of 

finite  length,  it  is  convenient  to  derive  the  theorems  of  the  correlation 

technique  for  continuous  records  and  infinite  length. 

The  auto-  and  crosscorrelation  functions  are  defined  as 


t     (x)  -  lim  ~~  J      f(t)f(t+t)dt  (B-l) 

t.-w   i     i 


and 


t. 
1    * 
♦.  (t)  =  lim   —  /   f(t)g(t+i)dt  (B-2) 

fe    t  -h.  2ti  _ti 


The  functions  f(t)  and  g(t)  are  the  records  under  consideration  and 
t,  the  argument  of  the  correlation  functions,  is  the  delay  time.   The 
autocorrelation  function  is  a  measure  of  the  similarity  of  a  record 
and  its  time  delayed  version,  while  the  crosscorrelation  function 
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relates  a  given  record  to  a  time  delayed  version  of  another  record. 
The  advantage  of  this  procedure  is  to  average  the  relationship  between 
the  records  over  the  length  of  the  records  and  thus  eliminate  any 
random  fluctuations  in  either  of  them. 

The  relationship  of  the  Fourier  transform  of  the  correlation 
functions  to  the  Fourier  transform  of  the  functions  is  of  considerable 
interest.   The  Fourier  transform  of  the  record  is 

F{f(t)}  -  f(M)  -  —  /"  f (t)elutdt  (B-3) 

/2ir  -°° 

The  Fourier  transform  is  a  complex  function  of  the  frequency   and 
let  it  be  noted  here  that 

f  (-11)  -  fCu)*  (B-4) 

where  the  asterisk  denotes  the  complex  conjugate. 

Let  the  Fourier  transform  representation  in  Eq.  (B-l)  be  substi- 
tuted for  the  function  f(t+x).   Reversing  the  order  of  integration 
and  noting  that 

lira  2F";_t  f(t)e"iUtdt  =  lim  5— I(-6>)  (B-5) 

t.-w   i    i  t.-"°   i 

1  x 

the  following  result  is  obtained 
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4»ff(T)   =  lira     -^- f2jMl(-bi)e   lt0Tdu  (B-6) 

t.-**3       i 


and  taking  the  Fourier  transform  on  both  sides 


J{»ff(T)}   =  lim^|f(u)]2  (B-7) 

t  .'-*»       i 


A  similar  result  is  obtained  for  the  cross  correlation  function, 
following  the  same  steps 


lim  ~  f  (-u)g(a,)  -  —  /*«-  (x)eluTdT  (B-8) 


From  these  expressions  the  possibility  of  evaluating  the  frequency 
response  of  a  system  is  obvious.   The  frequency  response  function  is 
defined  as 


«W-*^  (B-9) 


where  g(tii)  and  f(a)  are  the  Fourier  transforms  of  the  output  and  input, 
respectively.   Substituting  into  Eq.  (B-9)  according  to  Eqs.  (B-7)  and 
(B-8)  the  following  important  result  is  obtained 


,   .-,  v»2  F{i      (t)} 


iw     F{*ffiT)} 
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From  the  autocorrelation  functions,  only  information  on  the  amplitude 
of  the  frequency  response  function  is  available,  but  from  the  cross- 
and  autocorrelation  functions  the  complete  (complex)  frequency  response 
function  can  be  obtained, 

A  case  of  particular  interest  to  this  investigation  is  the  cross- 
correlation  of  a  harmonic  function  and  a  delayed  function  with  the  same 
period. 
Let       f(t)  =  A  sinwt 

g(t)  =  B  sin(wt-S)  (B-12) 


then 


**  to  "  ^T1  I   I  sinwt  sln(ui(t+T)-6)dt  (B-13) 

fg     2tf  -t, 


Here  the  crosscorrelation  of  a  limited  record  is  considered,  of  the 

length  2t .+t    .   Use  of  addition  theorems  for  the  harmonic  function 
l  max 

leads  to 


ff    to   "  irr~  f  ..[sin  cutCcosut   costS  +  sinux  sin5)  + 
8  i         1 


+  sinwt   coswt(sinojT    cos6   -   cosut   sinS)]dt  (B-14) 


Note  that 


i    2  1 

/  ,.  sin  lot  dt  =  t.  -  -r~   sin2ut. 

-t.  i   2u       i 
l 
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t. 

/   sinut  cosw.t  dt  =  0  (B-15) 

"  i 


and  the  result  is 


.   sin2(ut. 

*fg(T)  =  (I  "  tot  X)  AB  cos(a)t-'s)  <B_16) 


Thus  the  delaying  angle  6  is  preserved  in  the  crosscorrelation  function 
as  the  lag  of  the  cosine. 

For  the  autocorrelation  function  it  can  be  similarly  shown  that 


.    sin2ut.   „ 

fff(t)  =  (i — ^r^Acosat  <B-17> 


for       f  (t)  =  A  sin((ut+<5) 

From  the  crosscorrelation  and  the  autocorrelation  functions  both  the 
phase  lag  and  the  amplitude  ratio  can  be  determined,  as  shown  before 
for  a  more  general  case. 

It  should  be  noted  that  the  square  of  the  original  amplitude 
appears  in  the  autocorrelation  function.   If  the  original  record  con- 
tains several  harmonic  functions  the  ratios  of  their  amplitudes  in  the 
autocorrelation  function  is  the  square  of  the  amplitude  ratio  in  the 
record.   Thus  a  square  wave  with  the  period  —  ,  in  Fourier  series 
representation 


f  (t)  =  A(coswt  -  -j  cos2iot  +  t   cos3ut  -  .  .  .)  (B-18) 
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has  a  sawtooth  autocorrelation  function,  the  Fourier  representation 
of  which  is 


A2  1  1 

4   (t)  ■  —z   (coswt  +  —s   cos2ujt  +  —j  cos3un  +  .  .  .)     (B-19) 


Thus  the  harmonic  with  the  highest  amplitude  is  emphasized  in  the 
autocorrelation  function. 

The  advantage  of  the  use  of  the  correlation  function  instead  of 
the  direct  investigation  of  the  records  (by  frequency  analysis)  is 
that  the  noise  components  are  greatly  reduced  by  the  averaging  process 
of  the  correlation  function  computation,  which  requires  much  simpler 
operation  than  the  frequency  analysis.   This  is  illustrated  in  the 
following  analysis. 

Assume  a  record  containing  a  signal  and  a  noise  component 

f(t)  -  s(t)  +  n(t) 

the  noise  component  being  characterized  by  zero  expectation  of  the  mean 
and  certain  total  power  (variance) .   The  autocorrelation  function  in 
this  case  is 


1    'i 
i>ff(t)  =  lim  j~     /  *  (s(t)  +  n(t))(s(t+t)  +  n(t+x))dt   (B-20) 

t.-*»   i     i 


which  reduces  to  the  four  terms 
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i>   (t)  -  *   (t)  +  $      (t)  +  4   (t)  +  *  (t) 
f  f      ss      sn      ns      nn 


(B-21) 


If  the  noise  is  not  correlated  with  the  signal,  the  middle  two  terms 
approach  zero  and  only  the  signal  and  the  noise  autocorrelation 
functions  remain.   For  white  noise  and  infinite  integration  time  the 
noise  autocorrelation  function  is  a  {-function  in  i .      For  limited 
lengths  of  the  record  and  for  limited  bandwidth  the  autocorrelation 
function  of  the  noise  is  a  decaying  exponential.   In  the  case  of 
crosscorrelation  of  a  noisy  record  with  a  reference  record  that  does 
not  contain  noise  there  are  only  two  terms 


*„  (t)  =  *   (t)  +  *   Ct) 
fg      sg      ng 


(B-22) 


and  the  noise  component  is  eliminated  more  efficiently. 


B.2  Discrete  and  Finite  Correlation  Function. 

In  practical  applications  only  correlation  functions  of  finite 
length  are  obtained  and  most  often  in  a  discrete  (tabulated)  form. 
From  the  detailed  analysis  of  the  limitations  imposed  by  the  finite 
length  of  the  correlation  function  (7,  8)  two  important  conclusions 
shall  be  cited. 

The  finite  length  of  the  correlation  function  reduces  the  Fourier 
transform  to  the  Fourier  series  analysis,  with  only  discrete  frequencies 
present.   These  frequencies  are  at  a  distance  Af  =  —  if  the  total 
length  of  the  function  is  2T.  Two  frequencies  closer  than  Af  cannot 
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be  distinguished. 

The  discrete  nature  of  a  correlation  function  renders  it  impossible 
to  distinguish  between  a  frequency  component  and  another  component  of 
the  related  frequency  from  the  series 


f,  2kfN-f,  2kfN  +  f,  ...  ,  k  =  1,  2, 


where  f  =  1/2At  is  the  Nyquist  frequency.   To  avoid  the  ambiguity 
introduced  by  this  folding  process  it  is  preferable  that  the  amplitudes 
of  all  components  outside  the  principal  frequency  band,  0  to  f  ,  be 
negligible  compared  to  those  within. 
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APPENDIX  C:  Analysis  of  Correlator  Response 

C.l  Correlator  Amplitude  Response 

For  the  analysis  of  the  magnetic  tape  recordings  of  the  thermo- 
couple voltages  a  TMC  Correlation  Computer  (Model  257,  NE  1363)  in 
conjunction  with  the  TMC  Computer  of  Average  Transients  (CAT  400  C, 
NE  1364)  was  used.  In  the  preliminary  investigations  it  was  found 
that  the  amplitude  ratios  obtained  for  the  thermocouple  voltage  records 
did  not  match  the  phase  angle  differences  as  expected  from  the  theoret- 
ical model.   Since  the  phase  differences  agreed  with  the  expected  values, 
the  amplitude  responses  of  the  correlator,  the  tape  recorder,  the 
recording  and  the  playback  circuits  were  investigated.   The  correlator 
response  was  tested  by  three  methods:   autocorrelation  of  unbiased 
square  wave  input,  autocorrelation  of  unbiased  sinusoidal  input  and 
correlation  of  sinusoidal  inputs  with  various  bias  voltages,  periods 
and  correlator  settings. 

With  the  square  wave  input  17  runs  were  performed  with  flat-to- 
flat  amplitude  ranging  from  .162  V  to  5.57  V,  the  maximal  permissible 
input  voltage  of  the  correlator  being  +  3  V.   The  square  wave  was 
generated  by  an  oscillator  (Krohn-Hite  Oscillator,  Model  400  C,  NE  1343) 
and  the  flat-to-flat  voltage  was  determined  by  a  precision  voltmeter 
(Keithley  Electrometer  600  A,  NE  906) .   The  frequency  of  the  oscillator 
was  .148  Hz  and  the  correlator  settings  were:   A*A  (autocorrelation), 
128  points,  multiplier  16.   Timing  was  manual  and  an  integration  time 
of  2  minutes  was  used. 
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The  autocorrelation  function  obtained  was  fitted  by  the  analytic 
form 


y  =  ao  +  a1(-l)lnt(2x/a2)(l-(2x)moda2)  (C-l) 


where  a  is  the  constant  component,  a.,  the  amplitude  and  a?  the  period 
of  the  sawtooth  function.   The  sum  of  error  squares  was  minimized  by  a 
general  multivariable  function  minimization  code.   Subsequently  this 
code  was  found  to  be  too  slow  to  be  used  for  the  least  squares  analysis 
with  several  parameters  and  the  code  described  in  Appendix  F  was  developed 
for  the  analysis  of  the  experimental  data. 

The  results  of  the  square  wave  tests  are  plotted  in  Fig.  28.   On  a 
linear  plot  some  nonlinearity  was  visible  and  indeed  on  the  log-log  plot 
the  experimental  points  seem  to  define  a  line  with  exponent  (slope)  1.085. 
The  slope  was  determined  graphically  and  no  further  analysis  was  performed. 
as  the  square  wave  input  is  not  representative  of  the  records  analyzed 
on  the  correlator. 

A  second  set  of  correlator  amplitude  response  tests  was  performed 
using  a  sine  wave  input  generated  by  the  oscillator  (same  as  above).  The 
input  amplitude  was  again  measured  by  the  voltmeter.   15  tests  were 
performed  with  identical  inputs  and  8  tests  with  different  amplitudes 
to  the  two  correlator  inputs.   The  correlator  settings  were:  A*B 
(cross-correlation),  128  points,  multiplier  16.   The  period  of  the 
input  oscillations  was  10  seconds.   The  correlator  output  was  analyzed 
by  the  program  described  in  Appendix  F  and  the  results  are  plotted  in 
Fig.  29. 
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Flat-to-flat  voltage  (volts) 
Figure  28.   Correlator  amplitude  response.   Square  wave  input 
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Figure  29.   Correlator  amplitude  response  test.   Sine  »« 


ave  input. 
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Note  that  the  product  of  the  two  input  amplitudes  is  used  as  the 
abscissa  and  the  correlator  output  amplitude  as  the  ordinate,  in  Fig.  29. 
The  square  roots  of  these  quantities  are  related  for  the  first  correlator 
test. 

The  nonlinearity  is  again  evident.   To  determine  whether  it  is 
significant  a  regression  analysis  for  13  runs  with  equal  input  (auto- 
correlation runs)  was  performed.   Two  models  were  fitted  to  the  data  points, 
Model  A 

y  =  a  x2  (C-2) 

J  o 

Model  B 

y  =  a  xal  (C-3) 

*    o 

where  y  is  the  correlator  output  amplitude  and  x  is  the  input  voltage. 
A  simple  computer  program  Was  devised  for  the  sum  of  squares  minimization 
for  the  nonlinear  model  B.   For  model  A  the  least  sum  of  error  squares 
was  found  to  be  2.34  x  10   (arbitrary  units)  and  for  the  model  B  with 
exponent  2.172  the  sum  of  squares  was  3.34  x  10  .  The  99%  confidence 
margin  around  the  smaller  sum  of  squares  is  3.88  x  10  and  thus  the 
difference  in  exponent  can  be  considered  highly  significant.   At  the 
same  time  the  spread  of  the  points  around  this  curve  is  small  enough 
that  the  amplitude  information  could  be  used  for  the  measurements.   Yet 
the  eight  test  runs  with  different  inputs  were  considerably  off  the 
line  fitted  to  the  autocorrelation  points.   Invariably  the  points  for 
runs  where  correlator  input  B  had  higher  amplitude  than  A  are  above  the 
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line  and  vice  versa.   Note  the  very  good  agreement  between  the  exponent 
obtained  in  this  test  and  the  one  obtained  graphically  for  the  square 
wave  input  tests  (2.172  versus  2  x  1.085). 

To  establish  a  relationship  which  would  take  into  account  the 
effects  of  inputs  A  and  B  separately,  a  third  series  of  tests  was  per- 
formed. To  each  input  various  bias  voltages  were  applied  and  two 
periods  of  oscillations,  two  correlator  settings  and  two  integration 
times  were  used.   The  input  amplitudes  as  the  main  variables  of  investi- 
gation and  the  background  variables  xrere  randomized  for  the  30  test  runs 
performed.   The  oscillator  was  used  to  generate  the  sinusoidal  input 
voltage  and  the  exact  ratios  were  obtained  on  the  analog  computer,  the 
components  of  which  are  within  .1%  tolerance.   The  results  of  the  tests 
are  plotted  in  Fig.  30.   The  abscissa  is  the  product  of  the  two  input 
voltage  amplitudes  and  the  ordinate  is  the  correlator  output  amplitude. 
The  correlator  output  was  analyzed  with  the  computer  program  described 
in  Appendix  F. 

A  very  wide  spread  of  the  test  results  is  evident.  Although  some 
spread  was  expected  from  the  results  of  the  previous  test  for  different 
amplitudes  of  the  two  correlator  inputs,  a  detailed  inspection  reveals 
that  it  is  due  here  not  only  to  the  different  response  of  the  two  inputs 
but  also  to  a  considerable  extent  to  the  other  variable  conditions  in 
the  tests.   In  Fig.  31  the  10  test  results  are  plotted  for  which  the 
product  of  the  two  input  voltage  amplitudes  was  1.   The  spread  within 
the  group  of  four  tests  where  B  input  was  2  volts  is  more  than  10%, 
which  is  much  in  excess  of  the  tolerable  error.   The  bias  voltage  seems 
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Figure  30.   Correlator  amplitude  response.   Randomized  sine  inputs. 
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Legend  to  point  designation 

(V  V  a3>  V  a5) 
a  -  A  input  bias  (V) 

a  -  B  input  bias  (V) 

a^  -  oscillation  period  (sec) 

a^  -  correlator  delay  multiplier 

a^  -  integration  time  (min) 
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Figure  31.   Correlator  amplitude  response.   Sinusoidal  inputs,  product 
of  voltage  amplitudes  =1. 
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to  be  quite  important  though  not  the  only  influence  involved.   No  simple 
correlation  between  the  parameters  is  evident  and  even  if  a  satisfactory 
model  would  be  developed  it  would  be  very  difficult  to  apply  it  as  in 
the  experimental  temperature  records  the  bias  is  variable  for  both  inputs. 

At  this  point  it  was  concluded  that  the  amplitude  information 
obtained  from  the  correlations  of  the  experimental  records  cannot  be 
relied  upon  and  only  the  phase  information  was  used  in  the  experimental 
data  analysis. 

C.2  Response  of  Individual  Correlator  Channels. 

A  number  of  correlator  outputs  from  the  experimental  records  were 
analyzed  and  it  was  observed  that  the  randomness  criterion  for  the 
residues  (cf.  App.  F  2.7)  was  rarely  met.   The  residues  were  plotted 
and  a  recurring  pattern  was  observed,  similar  to  the  correlator  bias 
response,  Fig.  24,  obtained  later.   Since  the  prominent  features  of 
the  residuals  pattern  are  easily  identified  with  certain  correlator 
channel  numbers  the  effect  was  attributed  to  the  correlator.   Several 
tests  were  performed  with  a  constant  voltage  as  the  input  of  the 
correlator.   The  correlation  function  had  the  shape  as  presented  in 
Fig.  24,  and  not  only  was  the  general  trend  repeated  in  different  runs 
but  also  the  fine  detail  of  the  function,  i.e.  the  fluctuations  of  the 
output  for  individual  channels  is  not  random  but  due  to  some  character- 
istics of  the  equipment.   Further  it  was  found  that  the  amplitude  and 
the  sign  is  correlated  with  the  input  A.   It  is  guessed  that  the  cause 
of  this  unevenness  is  in  the  temporary  storage  for  the  input  A,  as  this 
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is  the  input  that  is  delayed  by  the  correlator.   The  amplitude  of 

the  largest  deviation  from  the  mean  horizontal  line  is  approximately 

-3 

.7  x  10   of  the  norm  count  for  a  1  V  bias  voltage  input. 

For  further  analyses  of  the  correlator  outputs  this  correlator 
background  function  was  incorporated  as  one  of  the  component  functions 
in  the  least  square  fit  in  the  code  described  in  Appendix  F.   The 
correlator  background  function  is  found  to  be  a  statistically  signifi- 
cant component  of  the  correlator  output  in  a  great  majority  of  cases. 
A  correlation  also  exists,  as  expected,  between  the  mean  deviation 
from  the  norm  count  and  the  correlator  background  function  coefficient. 

In  approximately  25%  of  the  cases  the  randomness  criterion  for 
the  residues  after  the  least  squares  analysis  (cf.  App.  F  2.7)  is  still 
not  satisfied.  A  possible  reason  for  this  may  be  that  the  background 
function  is  also  nonlinear ly  dependent  on  the  bias  of  the  input. 

C.3  Accuracy  of  Phase  and  Period  Determination, 

The  accuracy  of  the  phase  angle  determination  from  the  correlator 
output  can  be  estimated  from  the  phase  angle  found  in  the  correlator 
tests,  where  the  inputs  were  of  the  same  or  opposite  phase  with  no 
lag.   Of  the  30  tests  in  the  last  group  of  correlator  tests  (see 
preceding  section)  18  were  with  opposite  phase  and  12  with  phase  0. 
The  mean  deviation  from  the  expected  value  was  .00291  and  .00107, 
respectively  for  the  two  subgroups,  and  the  standard  deviation  estimates 
were  .00410  and  .00506.   It  can  thus  be  concluded  that  the  phase  angle 
error  shall  exceed  .01  radians  or  approximately  .5°  in  less  than  90% 
of  the  cases . 
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The  least  squares  analysis  of  the  correlator  output  sometimes 
gives  narrower  confidence  limits,  down  to  .001  rad,  but  this  only 
indicates  that  the  data  define  the  angle  with  that  precision  and  not 
that  these  are  the  limits  within  which  the  real  angle  is  located  with 
the  specified  probability. 

To  determine  the  accuracy  of  the  delay  increment  timing  of  the 
correlator  a  test  was  made  with  square  wave  input.  The  period  of 
oscillation  of  the  square  wave  was  measured  with  an  electronic  timer 
(Eckman  Time  Interval  Meter,  Model  7250R,  NE  763).   The  measured 
period  of  the  input  was  compared  to  the  period  of  the  autocorrelation 
function.   It  was  found  that  the  correlator  delay  timing  is  1.5  to  2% 
slow,  i.e.  the  delay  increments  are  that  much  longer  than  specified. 
In  present  measurements  the  period  of  the  temperature  oscillations 
was  determined  by  the  synchronous  motor  speed  and  the  gear  ratio  and 
not  from  the  correlator  output. 
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APPENDIX  D:  Analysis  of  the  Recording  Circuit,  Playback  Circuit  and 
Tape  Recorder  Response. 

D.l  Recording  and  Playback  Circuits. 

To  eliminate  high  frequency  electronic  noise,  filtering  circuits 
were  used  during  recording  of  the  thermocouple  voltages  and  during 
the  record  analysis.   The  signals  were  also  amplified  both  during 
recording  and  during  analysis.   For  the  filtering  and  amplifying 
circuits  the  components  of  an  analog  computer  were  used  (Electronics 
Associates  Co.   PACE  TR-10 ,  HE  816).   The  resistors  and  the  capacitors, 
which  determine  the  amplification  and  the  time  constants  are  guaranteed 
to  be  within  .1%  of  the  nominal  value.  Measurements  of  the  resistance 
of  the  resistors  showed  even  less  deviation  between  individual  parts. 
The  schematics  of  the  circuits  are  shown  in  Fig.  32. 

The  filtering  function  is  performed  by  an  integrator  with  resistance 
in  the  feedback  loop.   The  transfer  function  of  such  a  circuit  is 


U      R,/R. 

-^  =  f  1  (D-l) 

U.   Efcfs+1 


where  index  f  denotes  feedback  and  i  input  components.   The  direct 
current  gain  of  such  filter  is  Rf/R.  and  the  break  frequency  (angular) 
is  ub  =  RfCf. 

In  the  recording  circuit  two  filtering  units  were  used,  with  the 
following  eler.ents: 
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Figure  32. a.   Recording  filtering  and  amplifying  circuit. 
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Figure  32. b.   Playback  filtering  and  amplifying  circuit. 
Figure  32.   Filtering  and  amplifying  circuits. 
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C  -  10~5f 

Rf  =  io5  a 
r.  -  10  a 


Thus  the  filter  had  a  10  fold  amplification  and  the  time  constant 
was  1  sec,  the  corresponding  break  frequency  is  id,  =  1  sec   ,  f  =  .159  Hz. 
The  record  was  played  back  at  16  times  the  recording  speed  and  the  tape 


output  was  passed  through  another  filter  with  time  constant  .05  sec 

1   in    -1  f1 
B.  «  20  sec  ,  f„ 

b  £■ 


(break  frequency  ay  »  20  sec"  ,  f.  =  3.18Hz).   At  the  recording  speed 


this  would  correspond  to  a  break  frequency  a,    -   1.25  sec   ,  f,  =  .198  Hz. 
The  shortest  period  of  the  temperature  oscillations  used  was  1.92  minutes, 
which  corresponds  to  a  frequency  of  .0087  Hz  and  is  more  than  a  decade 
below  the  break  frequencies,  in  the  region  where  the  response  can  be 
considered  the  same  for  all  practical  reasons  as  for  direct  current. 

More  important  than  the  absolute  lag  caused  by  the  filtering  circuits 
is  the  difference  between  the  circuits.  A  test  was  performed  to  determine 
this.  A  sinusoidal  voltage  was  passed  through  the  four  parallel  recording 
circuits,  recorded  and  the  record  analyzed  on  the  correlator.   To  simulate 

the  actual  conditions,  different  amplifications  were  used  in  the  four 

3      5 
circuits,  ranging  from  10   to  10  ,  the  input  signal  being  previously 

appropriately  reduced  for  each  circuit.   The  correlator  outputs  were 

analyzed  by  the  computer  program  described  in  Appendix  F.   The  phase  angle 

determined  was  in  all  cases  within  .008  radians  (approx.  .45°)  from  0 

or  i,  depending  on  the  sign  of  the  record.   This  is  within  the  limits  of 

the  accuracy  of  the  phase  angle  determination  from  the  correlator  output 

(cf.  section  C  3.).   The  recording  and  playback  circuits,  and  the  tape  recorder 
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thus  do  not  seem  to  add  any  error  to  the  phase  between  the  Inputs. 
The  limits  of  accuracy  for  the  electronic  equipment  presently  used 
can  be  considered  .01  rad  or  .5°. 

D.2  Tape  Recorder  Amplitude  Response. 

In  the  examination  of  the  results  of  the  circuitry  and  the  tape 
recorder  phase  angle  accuracy  it  was  noticed  that  the  amplitudes  of 
the  correlation  functions  were  not  identical  as  they  should  have  been. 
As  a  possible  source  of  the  incongruence  the  tape  recorder  amplitude 
response  was  investigated.  A  sinusoidal  voltage  was  introduced  to  the 
inputs  of  the  four  recording  channels.   The  autocorrelation  function 
of  each  record  was  obtained  and  subsequently  analyzed  with  the  computer 
program.  The  results  are  presented  below. 

Tape  recorder       Normalized  autocorrelation      Square  root  of 
channel  function  amplitude        normalized  amplitude 

1  .03550  .188 

2  .04600  .214 

3  .04640  .215 

4  .0496  .222 

The  differences  are  due  to  the  imperfect  condition  of  the  recorder 
data  converter  circuits.   The  weaker  gain  of  channel  1  was  observed 
early  in  the  course  of  investigations,  but  the  repair  was  not  possible  as 
the  malfunctioning  component  substitute  was  not  available.   The  imperfec- 
tion of  the  tape  recorder  further  justifies  the  decision  not  to  attempt 
to  use  the  amplitude  information  of  the  temperature  oscillations. 
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APPENDIX  E:   Regression  Analysis 

E.l  Least  Squares  Principle 

Of  several  possible  criteria  for  obtaining  a  "best"  fit  of  an 
analytic  curve  to  a  set  of  data  points  the  least  squares  principle  is 
most  often  applied.   In  the  case  of  normally  distributed  errors  the 
curve  fitted  in  this  sense  is  the  curve  of  maximal  likehood,  as  can 
be  seen  from  the  following  (14) . 

Assume  that  the  errors  of  a  sample  y  ,  y_,  .  .  . ,y  ,  corresponding 
to  the  independent  variable  vectors  X  ,  X_,  ...  X  from  the  time 

values  are  random  with  expectation  zero,  with  no  correlation  (covariance 

2  2 

matrix  is  diagonal  V(e  =  la  )  and  normally  distributed,  e-N(0,Ia  ). 

The  likehood  function  for  the  sample  is  in  that  case 


P(yr  y2,  ....  y  )  -  H     \/2  ^^V2^   = 

i=l  0"(2n) 


1     exp(-eTe/2o2)  (E-l) 


o  (2?r) 


The  maximum  likehood  estimate  of  an  analytical  model  is  thus  the 
one  that  minimizes  the  sum  of  squares  of  the  differences  between  the 
observed  values  and  the  predicted  values.   If  there  is  reason  to  believe 
that  the  errors  have  a  different  distribution  than  normal  other  criteria 
for  the  maximal  likehood  estimate  can  be  obtained. 
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E.2  Least  Square  Determination  of  Linear  Parameters 

A  set  of  observations  is  given,  y.  ,  y„,  .  .  . ,  y 
and  the  assumed  analytical  model  is  a  linear  combination  of  k  functions 
of  the  independent  variables 


y(x)  =  bxZ(x)  +  b2Z(X)  +  .  .  .  +  bkZ(X) 


(E-3) 


The  problem  of  finding  the  values  of  the  linear  constants  b  is  one  of 
solving  the  system  of  equations,  in  matrix  notation 


w   w W" 

W   Z2(X2)    W 


W    W 


z,  (x  )    b,     y 

k  n       k      '■" 


(E-4) 


Zb  =  y 


(E-4a) 


In  the  practical  cases  when  the  observations  y.  are  subject  to 
errors  the  number  of  the  observations  should  largely  exceed  the  number 
of  the  parameters  to  be  determined  and  we  have  an  overdetermined  incon- 
sistent system.  A  consistent  system  is  obtained  by  adding  an  error 
vector, 


Zb  =  y  +  e 


(E-5) 


and  the  solution  of  this  overdetermined  system  shall  be  sought  in  the 
least  sum  of  square  errors  sense. 

The  sum  of  squares  of  the  errors  is 


eTe  =  (Zb-y)T(zb-y)  =  bTzTZb  -  yT  zb  -  bTzTy  +  yTy 


The  necessary  condition  for  a  minimum  is  that  all  the  derivatives 
according  to  individual  parameters  are  equal  to  zero 


T 
3(^be)  -  t0,  0,  ...,0,  1,  0,  ...,0jZTZb  +  bTZTZ,0,  ...,0,  1,  0, 
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the 

(E-6) 

,.,0jT  -  2^0,  ...,  0,  1,  0,  ....  0jZTy  =  0     (E-7) 


The  k  such  equations  for  the  parameter  vector  b  at  the  point  of  minimal 
sum  of  squares  can  be  assembled  into  a  new  matrix  equation 


ZTZb  =  ZTy  (E-8) 


T 
If  the  square  matrix  Z  Z  is  not  singular  we  have  the  solution 


bQ  -  (ZTZ)1  ZTy  (E-9) 


which  satisfies  the  system  Eq.  (E-4)  in  the  sense  of  minimizing  the 
sum  of  squares  of  the  errors.  At  this  point  the  standard  deviation 
of  the  errors  can  be  estimated,  the  unbiased  estimate  being 
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s2  =  eTe/(n-k)  (E-10) 


To  check  the  validity  of  the  assumed  analytical  model  the  errors 
(residues)  should  be  tested.   If  the  errors  are  not  random  the  analy- 
tical model  Eq.  (E-3)  is  not  appropriate. 

E.3  Confidence  Limits  of  the  Parameters 

If  the  errors  are  distributed  normally  the  F-test  is  applicable 
to  the  sum  of  squares  of  the  errors.   The  minimal  sum  of  squares  of 
the  differences  between  the  observations  and  the  analytical  model  points 
is  compared  to  the  sum  of  squares  of  the  differences  between  the  least 
square  analytical  model  (parameter  vector  b  )  and  the  analytical  model 
with  a  different  parameter  vector,  b. 

Here  the  error  sum  of  squares  of  the  observations  around  the  least 
square  solution  has  n-k  degrees  of  freedom  and  the  additional  sum  of 
squares  of  a  compared  fit  has  k  degrees  of  freedom  (n  being  the  number 
of  data  points  and  k  the  number  of  parameters  in  the  model  equation) . 

The  two  sums  of  squares  are  defined  as 


ss(b|bo)  =  (y[b]  -  y[bQ])T(y[b]  -  y[bQ])  =  (b-bo)^TZ(b-bQ) 


(E-ll) 


SS(b  )  =  (y-y[bn])T(y-y[bn])  -   (y-Zb  )T(y-Zb)  (E-12) 


where  y  are  the  observed  values  and  y  (b  )  and  y.(b)  are  the  values 
predicted  by  the  model  with  different  parameter  vectors.   The  criterion 
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for  the  nonsignificance  of  the  additional  error  sum  of  squares  at  the 
(l-o)  confidence  level  is 


SS(b]b  )/SS(b  )  <  k  F,n   ,(k,n-k)/(n-k)  (E-13) 

1  o      o  -    {!■-&) 


This  relation  can  be  expressed  as 


(b-b  )TZTZ(b-b  )  <  (Zb  -  y)T(Zb  -y)  k  Fn  .  (k,n-k)/n-k)  =  ASS 
o        o  -    o         o       {±-aj 

(E-14) 


In  the  k-  dimensional  space  of  the  parameters  b .  the  relation  defines 
a  region  inside  which  the  actual  parameter  vector  is  situated  with  the 
specified  probability,  on  the  basis  of  the  information  contained  in  the 
data.   The  equality  sign  defines  the  confidence  region  limit.   As  we 
have  a  quadratic  form  on  the  left  side  the  form  of  the  confidence 
region  is  a  hyperelipsoid.   For  visual  presentation  the  case  of  two 
linear  parameters  shall  be  used  (Fig.  33). 

The  confidence  region  is  fully  determined  by  the  vector  b,  the 

T 
matrix  Z  Z  and  the  confidence  margin  of  the  sum  of  squares,  ASS.   For 

each  individual  parameter  a  confidence  interval  can  be  obtained.   The 
confidence  interval  for  each  parameter  can  be  defined  either  as  the 
width  of  the  elipsoid  on  the  line  parallel  to  the  axis  of  the  parameter 
(points  C  and  D  in  Fig.  33),  or  by  the  "support  plane"  distance,  i.e. 
by  tangential  planes  on  the  elipsoid  (points  A  and  B) .   The  later  defin- 
ition is  more  suitable  as  such  interval  contains  the  whole  range  that 
the  parameter  can  assume  within  the  confidence  region.   The  two  intervals 
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coincide  when  the  principal  axes  of  the  elipsoid  are  parallel  to  the 

T 
axes  of  the  parameters,  i.e.  when  the  matrix  Z  Z  is  diagonal  (the  functions 

Z.(X)  are  orthogonal).   Note  that  the  support  plane  confidence  limits 

should  not  be  applied  simultaneously  to  several  parameters,  e.g.  point 

E  in  Fig.  33  lies  outside  the  confidence  region  though  both  parameters 

are  within  the  confidence  intervals. 

The  extreme  value  of  a  parameter  on  the  confidence  region  boundary 

is  reached  in  the  point  that  satisfies  Eqs.  (E-15)  and  (E-16) .   The 

equation  of  the  confidence  region  boundary  is 


(b-b  )TZTZ(b-b  )  -  ass  (E-15) 

o        o 


and  the  partials  according  to  other  parameters  are 


a/3b.[(b-b  )TZTZ(b-b  )]  =  0        i  =  1,  2,  ...,  k  lft 

1      o  o 

(E-16) 

After  differentiation  k-1  equations  are  obtained 


t0,  0,  ....  0,  1,  0 Oj  ZTZ(bmj-bQ)  -  0  CB-17) 


where  the  only  element  different  from  zero  in  the  first  vector  is  at 
position  j . 

Substituting  these  equations  into  Eq.  (E-15)  we  have 

,0,  ...,  0,  b.   .,  -  b.,  0,  ...,  0,  ZTZ  (b  ,-b  )  =  ASS    (E-18) 

1  j  ,mj     j  i  mj   o 
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We  can  treat  the  equations  for  the  support  plane  points  of  all  parameters 
together  by  forming  a  matrix  equation 


bl,mrbl,0 


0 

0 

b2,m2-b2,0        • 

0 

.    . 

ZTZ 

b     -b 
ml     o 

b .  -b 

mk     o 

0 

■    ■        bk,mk"bk,0 

ASS  '  I 


(E-19) 


where  b   indicates  the  parameter  vector  at  the  support  plane  point  for 
parameter  i  and  b.   .  indicates  the  value  of  parameter  i  at  such  point. 

r  1,511 

Multiplying  by  corresponding  inverses  and  noting  that  the  multipli- 
cation with  diagonal  matrices  is  commutative  we  obtain 


b  .-b 

ml   o 


b  ,-b 

mz  o 


b  ,-b 

mk  o 


b-i   i~bi  n   ° 
l,ml  1,0 

2,m2  2,0 


bk,mk"bk,0 


(Z^)1  ASS 


CE-20) 


We  are  at  this  point  interested  only  in  the  confidence  limits  of  the  - 
individual  parameters  and  not  in  the  parameter  vectors  in  the  support 
plane  contact  points.   For  the  confidence  limits  of  the  individual 
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parameters  we  have 


(b.   ,-b.  „)  »  c.  .ASS 
1,1111  1,0      1,1 


(E-21) 


t-,.1 


where  c    is  the  diagonal  element  of  the  matrix  (Z  Z)  .  The  half 


width  of  the  support  plane  confidence  interval  is 


b4      .-b.  .  -  [c.  .eTe  k  Fn   ,(k,n-k)/(n-k)]1/2  (E-22) 

i,mi  1,0     1,1       (1-a) 


E.A  Linearized  Confidence  Limits  for  Nonlinear  Parameters 

The  model  which  is  fitted  to  the  observed  data  may  contain  nonlinear 
parameters.   The  matrix  method  described  above  is  not  applicable  for 
finding  the  least  squares  solution  nor  is  the  method  for  obtaining  the 
confidence  limits.   However  the  model  function  can  be  expanded  in  a 
Taylor  series  and  only  the  linear  term  retained  if  the  confidence  region 
is  small  enough  that  this  approximation  is  sufficiently  accurate.   The 
additional  error  due  to  the  deviation  in  the  parameters  is  then  given  by 


3bx  tXlJ     3b2  "V 


dbx   °V    3b2  lV 


fh   <XJ 

9b-   n 


3b    n 


k 


^-   (x  ) 
3b,  lV 
k 


&-  (x  ) 

3  b.  UnJ 
k 


(b-b  )  =  P(b-b  ) 


(E-23) 


124 


T  T 

and  the  additional  sum  of  squares  is  (b-b  )  P  P(b-b  ).   The  approximate 

confidence  region  is  limited  by  the  surface 


(b-b  )TPTP(b-b  )  =  ass 


(E-24) 


A  similar  procedure  can  be  followed  in  determining  the  support  plane 
confidence  interval  as  in  the  linear  case  and  the  result 


(b.   .-b,  J   =  Y.  .ASS 
J.mj   j,0     'j,j 


(E-25) 


T   I 
is  obtained,  where  y .  .  is  the  diagonal  element  of  the  matrix  (P  P)  . 
J.J 


2,m2+ 


2,m2 


Figure  33 .   Parameter  confidence  region  for  two  linear  parameters . 
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APPENDIX  F.   Correlator  Output  Analysis  Computer  Program 

F.l  General 

From  the  correlator  output  the  prominent  harmonic  component  is 
to  be  obtained,  as  well  as  its  amplitude,  phase,  period  and  possibly 
higher  harmonic  components  with  periods  that  are  multiples  of  the 
fundamental  harmonic  period.   Because  the  records  also  contain  some 
drift,  usually  nonlinear,  and  furthermore  since  the  response  of 
individual  channels  of  the  correlator  is  not  uniform  (cf.  Appendix  C) , 
the  appropriate  analytic  form  which  has  to  be  fitted  to  the  correlator 
output  is 

npoly-1 
y,  -  K  +  >     b.PVx.)  +  b   ,  F  ,  (x.)  + 

7i    -1   1=1     3  J   i     npoly  cb  l 


nharm  .  .  0  . 

V  r.  2TT-1X1   ,  ,  -    27TJX-:   ■,  /^  i  \ 

>  [b   ,  , .cos  — -J-*-  +  b         sin  — J-3-  ]  (F-l) 

'     npoly+j      tQ     npoly+j+1      tQ 


where  P'.  (x.)  is  a  polynomial  of  order  j  and  F  ,  (x.)  is  the  correlator 

bias  response  function.   The  Legendre  polynomials,  approximately 

orthogonal  on  the  sample  points,  are  used,  though  this  does  not  simplify 

the  computation  as  the  least  squares  fit  is  obtained  by  the  matrix 

method. 

One  parameter  in  the  expansion  appears  nonlinear ly,  the  period  of 

the  fundamental  harmonic,  t  .   If  several  nonlinear  parameters  were 

o  * 

involved  a  general  sum  of  squares  minimization  algorithm  would  have  to 
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be  used  (14,  23),  but  with  only  one  nonlinear  parameter,  a  simpler 
and  faster  method  can  be  applied. 

The  input  data  is  considered  to  be  an  equidistant  table  of 
functional  values,  and  only  the  ordinate  at  each  point  is  supplied. 
The  abscissae  are  generated  and  associated  internally  with  each  data 
point.   This  is  necessary  because  erroneous  points  are  rejected  by 
the  code  and  the  data  set  is  no  longer  equidistant. 

The  code  features  the  ability  to  determine  the  initial  guess  for 
the  prominent  frequency  period.   The  initial  guess  has  to  be  close 
enough  to  the  true  value  as  otherwise  a  local  minimum  of  sum  of  squares 
may  be  obtained  instead  of  the  desired  solution.  The  recognition  of 
the  prominent  harmonic  is  possible  even  if  several  erratic  points  are 
present  or  if  due  to  general  slope  the  local  extremes  of  the  harmonic 
component  are  inconspicuous. 

Erratic  points  in  the  correlator  output  are  often  observed,  due 
to  malfunction  of  the  electronic  components  of  the  correlator  or  in 
the  transfer  of  data.   The  residuals  of  the  least  squares  fit  are 
checked  and  each  of  them  compared  to  the  standard  deviation  estimate. 
If  the  sum  of  squares  of  the  residuals  above  a  present  value  (3.291a 
for  99.9%  confidence  level)  exceeds  half  of  the  confidence  limit 
margin  such  points  are  rejected  and  the  fitting  of  the  analytical 
model  repeated. 

The  residuals  are  subsequently  checked  for  randomness  and  the 
result  of  the  check  is  displayed. 

Ultimately  the  confidence  limits  for  the  parameters  are  calculated. 
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For  the  nonlinear  parameter  only  the  linearized  confidence  limit 
estimate  is  obtained.   For  the  harmonics  the  amplitude  and  the  phase 
angle  is  calculated  as  well  as  the  confidence  limits  for  the  phase 
angle.   The  parameters  correlation  matrix  is  also  displayed. 

Several  additional  output  options  are  available,  various  matrices 
can  be  displayed  during  the  least  squares  analysis  and  the  residuals 
can  be  punched  on  cards  for  further  analysis. 

Details  of  the  algorithms  employed  are  given  below.   The  program 
is  written  in  FORTRAN  IV  and  double  precision  has  been  used  throughout. 

F.2  Code  Description 

The  program  provides  for  the  maximum  of  15  parameters.   If  the 
number  of  linear  parameters  alone  is  15  the  confidence  limits  estimate 
for  the  nonlinear  parameter  is  not  available.   The  present  coding 
provides  for  a  maximum  of  256  input  data  points . 

The  code  is  composed  of  the  following  programs 

Name  No.  of  FORTRAN  statements 

MAIN  232 

LINLSQ  46 

ZXIIN  (ZXI)  35 

ANGLE  8 

XIMIN  90 

PFRE  77 

RCKECK  32 

REJECT  36 
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Name              No.  of  FORTRAN  statements 

WMAT  55 

MATINV  76 

Total  687 

FORTRAN  library  subroutine  TIME  is  called  for  timing  of  the  execution 
of  each  case.  The  program  was  tested  and  operated  on  the  KSU  IBM  360/50 
system. 

F.2.1  Main  Program 

The  arrays  used  in  the  program  are  dimensioned  to  allow  up  to  256 
points  and  up  to  15  parameters. 

The  batch  and  case  data  are  read  as  specified  in  the  section  F.3. 
According  to  the  specified  value  of  the  variable  NPREL  the  REJECT  and 
PFRE  subroutines  are  executed  or  skipped.   If  NPREL  is  larger  than  one 
all  portions  of  the  main  program  are  executed,  as  described  in  the 
following. 

Segments  of  the  input  data  table  of  length  KPTRE  are  submitted 
to  the  subroutine  REJECT  for  inspection  for  erratic  points.   For  the 
first  quarter  of  the  first  segment  the  rejection  normal  deviate  is 
set  at  1.96  (95%  confidence).   In  subsequent  steps  only  the  rejection 
of  the  middle  two  quarters  of  the  segment  is  effective,  at  the  99% 
confidence  level,  and  the  segments  overlap  to  cover  the  whole  length 
of  the  data  series.   In  the  last  quarter  of  the  last  segment  the 
rejection  is  again  effective  at  95%  confidence  level. 

Subsequently,  if  NPREL  is  specified  greater  than  1  a  least  square 
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polynomial  of  the  order  NPKEL-1  is  subtracted  from  the  data  series. 
Matrix  method  for  the  calculation  of  the  least  squares  polynomial  is 
used.   This  option  is  provided  to  enable  the  recognition  of  the 
prominent  harmonic  even  if  it  is  superimposed  on  a  steep  slope.   In 
such  case  there  may  be  no  local  maxima,  on  the  basis  of  which  the 
subroutine  PFRE  calculates  the  initial  guess.  After  the  subroutine 
PFRE  is  invoked  the  original  data  series  is  restored  for  further 
processing. 

If  NDUMP  is  specified  larger  than  one  only  points  with  index 
divisible  by  NDUMP  are  retained.  This  truncation  of  data  enables 
faster  calculation  and  often  the  precision  lost  is  not  crucial. 

The  polynomial  columns  of  the  matrix  Z,  their  number  specified 
by  NPOLY  +  1,  are  initialized  for  the  points  that  were  not  rejected. 
The  independent  variable  (lag  time)  is  transformed  so  that  the  data 
points  are  on  the  interval  -1  to  1.   Column  2  of  the  matrix  Z  contains 
the  value  (i-n/2)/(n/2)  where  i  is  the  index  of  the  data  point  and  n 
is  the  number  of  the  data  points.   This  column  provides  the  argument 
for  the  calculation  of  the  harmonic  function  in  the  subroutine  ZXIIN 
(ZXI) .   The  other  polynomial  columns  of  the  matrix  Z  are  initialized 
with  Legendre  polynomials,  which  are  approximately  orthogonal  on  the 
discrete  set  of  points  in  the  interval  -1  to  1.   If  IDEBUG  is  divisible 
by  29  (cf.  F-3.1)  the  column  NPOLY+1  is  initialized  with  a  special 
function  which  has  to  be  provided  as  input  data  for  the  batch  of  cases. 
This  provision  enables  inclusion  of  the  correlator  bias  response 
function  as  one  of  the  components  of  the  correlation  function  expansion. 
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The  harmonic  functions  columns  are  initialized  by  the  subroutine 
ZXIIN  (entry  2X1). 

The  minimization  of  the  sum  of  squares  by  adjusting  the  period 
of  the  fundamental  harmonic  is  performed  by  the  subroutine  XIMIN  and 
the  parameter  vector  and  the  vector  of  the  residuals  is  returned. 
The  check  for  randomness  of  the  residuals  is  performed  by  the  subroutine 
RCHECK. 

For  the  determination  of  the  confidence  limits  of  the  parameters 
the  matrix  of  the  partial  derivatives  of  the  fitted  function  according 
to  the  parameters  is  formed  (cf.  section  E.4).   The  partials  according 
to  the  linear  parameters  are  the  values  of  the  component  functions  and 
only  the  derivative  according  to  the  nonlinear  parameter,  the  period 
of  the  fundamental  harmonic,  is  to  be  calculated.   If  the  number  of 
linear  parameters  is  equal  to  the  maximal  dimension  of  the  array  Z 
(presently  15)  the  confidence  limits  for  the  nonlinear  parameter  are 
not  available.   The  confidence  limits  are  determined  according  to 
Eq.  (E-22) .   The  Fpt.„(k,  120)  values  are  supplied  as  the  code  data, 
where  k  is  the  number  of  parameters,  1  to  15. 

The  residuals  are  checked  for  erroneous  points.   The  rejection 
level  is  set  at  99.8%  confidence  (normal  deviate  3.09)  and  the  criterion 
for  repeating  the  least  square  analysis  is  the  sum  of  squares  of  the 
residuals  above  normal  deviate  3.291.   If  this  sum  exceeds  half  of 
the  sum  of  squares  confidence  margin  at  95%  confidence  level  (Eq.  (E-14)) 
the  calculations  are  repeated. 

The  amplitude  of  each  harmonic  is  calculated  from  the  cosine  and 
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sine  expansion  parameters,  as  well  as  the  phase  lag  in  radians  and 

independent  variable  increments  (lag  units),  so  that  the  harmonic  can 

2tt 
be  represented  as  A  cos  (—  x~<j>),  A  being  the  amplitude  and  <f  the 

o 

phase  lag.   If  the  confidence  interval  of  both  components  of  the 

harmonic  (cosine  and  sine)  includes  the  origin  the  amplitude  is 

considered  insignificant  and  is  not  calculated.  A  pessimistic  estimate 

of  the  phase  angle  confidence  limits  is  obtained  by  taking  the  larger 

of  the  confidence  intervals  of  the  two  parameters  at  the  distance  of 

the  amplitude  of  the  harmonic  as  the  angular  confidence  interval. 

T   I 
The  parameters  correlation  matrix  is  the  matrix  (P  P)  ,  where  P 

is  the  matrix  of  the  partial  derivatives,  Eq.  (E-2  3),  normalized  so 

that  each  row  and  column  is  divided  by  the  square  root  of  the  diagonal 

element.   If  the  component  functions  are  mutually  orthogonal  on  the 

given  set  of  points  the  parameter  correlation  matrix  is  a  unit  matrix. 

Subroutines  needed  by  this  program  are  listed  above  and  their 
description  is  given  in  the  following  sections. 

The  output  format  is  illustrated  by  two  sample  outputs  presented 
in  sections  F.6  and  F.7.   The  effect  of  the  IDEBUG  options  on  the 
output  is  described  in  section  F.3.1. 

F.2.2  Subroutine  REJECT 

Given  a  vector  of  equidistant  functional  values  of  the  length 
M  (maximum  50)  the  subroutine  REJECT  constructs  by  the  matrix  method 
an  approximation  polynomial  of  the  specified  order  N-l  (N  maximum  is 
10) .   The  deviations  of  the  input  data  points  are  compared  to  the 
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standard  deviation  estimate  for  the  deviations  around  the  approximation 
polynomial.   If  the  ratio  exceeds  the  specified  value  (FSS)  the  data 
point  is  set  equal  to  zero,  otherwise  it  is  set  equal  to  one. 

If  in  the  subsequent  references  of  the  subroutine  the  length  of 
the  functional  value  vector  and  the  order  of  the  polynomial  are  the 
same  the  matrix  operations  of  the  least  squares  fitting  are  not  repeated. 

Subroutine  MATINV  and  subroutine  package  WMAT  are  referred  by 
this  subroutine. 

F.2.3  Subroutine  PFRE 

Given  a  vector  of  equidistant  fur.ctional  values  the  amplitude, 
phase  and  period  of  the  prominent  frequency  is  estimated.   Local 
maxima  and  minima  of  the  function  are  determined,  the  criterion  being 
that  a  given  point  is  higher  (lower)  than  any  of  the  four  neighbouring 
points  on  either  side.   The  most  frequent  distance  between  adjacent 
maxima  and  minima  is  taken  as  the  estimate  of  the  period.   Phase  is 
determined  from  the  location  of  the  extremal  points  and  the  amplitude 
estimated  as  one-half  of  the  difference  between  the  average  of  the 
maxima  and  minima.   In  order  that  the  subroutine  be  effective  the 
maxima  must  be  recognizable  and  any  erroneous  points,  which  may  be 
mistaken  for  extremes  must  be  identified.   The  functional  value  vector 
must  be  previously  scanned  by  the  REJECT  subroutine  and  if  overall 
slope  is  excessive  a  polynomial  has  to  be  subtracted. 

F.2.4   Subroutine  XLMIN 

In  this  subroutine  the  minimization  of  the  sum  of  squares  for  the 
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nonlinear  parameter,  the  period,  is  performed.  A  good  approximation 
for  the  period  has  to  be  given,  sufficiently  close  to  the-  minimizing 
value.   Two  points  are  chosen,  one  on  each  side  of  the  guess  at  the 
distance  6  ,  the  convergence  criterion  for  the  period.   Sums  of  squares 
are  calculated  for  all  three  points  by  the  subroutine  LINLSQ,  as  well 
as  the  least  squares  parameters  vector  and  the  residuals.   The  dependence 
of  the  sum  of  squares  on  the  period  is  approximated  by  a  parabola  through 
the  three  points.   The  point  at  the  bottom  of  the  parabola  is  taken  as 
a  new  guess  and  the  point  with  the  maximal  sum  of  squares  is  rejected. 
The  first  two  criteria  for  convergence  are: 

!•   (x.~X-)  <  1j  X  being  the  deviation  of  the  period  estimate  from  the 
initial  guess  in  the  units  of  5  ,  the  convergence  criterion 
specified  on  the  case  input  or  assigned  by  default  (.05  delay 
increments)  in  the  calling  program. 
2.   (SS(Xi)  -  ss(X.))max  <  <5fASS  kF(1_a)(k,n-k)/(n-k)  i.e.,  the  maximal 
difference  of  the  sums  of  squares  must  be  a  fraction  of  the  confi- 
dence limit  margin.   6__  is  specified  in  the  case  input  or  assigned 
by  default  (6  =  .05)  by  the  calling  program. 
When  these  convergence  criteria  have  been  reached  it  is  further 
required  that  the  middle  point  in  x  be  the  one  with  the  lowest  sum  of 
squares.   The  three  point  set  is  moved  till  this  is  achieved.   Thus 
satisfactory  assurance  is  given  that  the  minimum  point  has  been  enclosed. 
By  specifying  19  as  a  factor  in  IDEBUG  this  last  convergence  criterion 
is  omitted. 

The  optional  printouts  from  this  subroutine  are  listed  in  F.3.1. 
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Following  subroutines  are  called  by  this  subroutine:   ZXIIH  (entry  ZXI) , 
L1NLSQ,  MATINV  and  the  subroutine  package  HMAT. 

F.2.5  Subroutine  ZXIIN,  ZXI 

The  columns  of  the  matrix  Z  containing  the  trigonometric  functions 
are  adjusted  to  the  current  period  of  the  fundamental  harmonic.   To 
save  computing  time,  sine  and  cosine  values  are  computed  directly  for 
each  10   point  on  the  equidistant  scale  and  for  the  intermediate 
points  the  values  are  calculated  by  the  addition  theorem  formulas.   As 
the  matrix  Z  does  not  contain  all  the  equidistant  points  the  value  in 
the  column  2  of  Z  is  used  as  the  argument  of  the  harmonic  function 
(see  section  F.2.1  on  the  values  in  polynomial  columns  of  matrix  Z). 

The  subroutine  is  dimensioned  for  the  maximal  size  of  the  matrix, 
Z  is  256*15. 

F.2.6  Subroutine  LINLSQ 

Given  a  matrix  of  the  values  of  the  component  functions  and  the 
observed  values  to  be  fitted  by  a  linear  combination  of  the  component 
functions,  the  vector  of  linear  parameters  (coefficients  of  the  component 
functions),  the  error  vector  (residuals)  and  the  sum  of  squared  errors 
are  calculated.   The  matrix  method  is  used,  as  described  in  Appendix  E. 

The  subroutine  arrays  are  presently  dimensioned  for  up  to  15 
parameters  and  up  to  256  observation  points,  i.e.  the  maximal  dimensions 
of  the  Z  matrix  are  256*15.   According  to  specified  IDEBUG  values  at 
each  execution  of  the  subroutine  the  arrays  2,1  1,(1  I)   and  the  error 
vector  can  be  displayed. 
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1.2.7  Subroutine  RCHECK 

The  randomness  of  the  residuals  is  tested  by  a  "run"  test;  the 
number  of  uninterrupted  sequences  of  residuals  of  equal  sign  is  compared 
to  the  expected  number  of  such  sequences  if  either  sign  would  be  equally 
probable  for  each  residual. 

If  r  is  the  number  of  runs  and  n  is  the  number  of  residuals  we 
have  (20) : 
Expected  value 

E(r)  =  (2n-l)/3  (F-2) 

Variance 

V(r)  =  (16ii-29) /30  (F-3) 

The  distribution  approaches  the  normal  distribution  for  large  n.  Normal 
distribution  is  assumed  and  if  the  normal  deviate  is  less  than  1.645  the 
randomness  is  not  rejected  at  the  95%  confidence  level. 

F.2.8  Subroutine  package  KMT 

Matrix  operations  are  performed  at  various  entry  points  of  the 
subroutine.   Variable  dimensions  of  the  arrays  can  be  used  and  the 
dimensioning  is  compatible  with  the  WATFOR  compiler. 

F.2.9   Subroutine  MATINV 

Matrix  inversion,  solution  of  a  system  of  linear  equations  and 
the  calculation  of  the  determinant  value  are  performed  by  the  subroutine. 
The  variable  dimension  feature  is  incorporated.   The  subroutine  was 
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written  at  Argorme  National  Laboratory. 

F.2.10  Function  ANGLE 

The  angle  in  radians  in  the  range  -.436  to  2ti-.A36  is  returned, 
given  the  non-normalized  cosine  and  sine  components  as  the  arguments. 

F.3.   Input  Preparation 

The  data  deck  consists  of  the  initialization  data,  case  data  sets 
and  end-of-batch  data. 

F.3.1  Initialization  Data 

Item  1.1     Format  2013;  FORTRAN  label  IDBG(I),  I  ■  1,20 

Product  of  IDBG(I)  is  formed  (default  for  each  IDBG(I)  ■  1) 
and  tested  in  the  code  for  the  content  of  a  particular 
prime  factor.   The  prime  factors  have  the  following  effect: 

2  The  polynomial  part  of  matrix  Z  is  printed  in  the  main 
program  (once) . 

3  Matrix  of  derivatives,  P,  is  printed  by  the  main  program 
(once) . 

5  Matrix  Z  is  printed  at  each  entry  in  subroutine  LINLSQ. 

T 
7  Matrix  Z  Zis  printed  at  each  entry  in  subroutine  LINLSQ. 

T  I  T-7 
11  Product  (1   Z)  (Z  Z)  is  printed  at  each  entry  in  subroutine 

LINLSQ  (check  for  inversion  accuracy). 
13  Error  (residuals)  vector  is  printed  at  each  entry  in 

subroutine  LINLSQ. 
17  Values  of  XINEW  and  FINEW  in  subroutine  XIMIN  are  printed 
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at  each  computation, 

19  Third  convergence  criterion  is  omitted  for  exit  from 
XIMIN  subroutine  (cf.  section  F.2.4). 

29  A  special  function  is  introduced  to  the  matrix  Z  in  the 
column  NPOLY  +  1,  the  function  thus  becomes  one  of  the 
component  functions  for  expansion  of  the  data  set.   If 
this  option  is  specified  the  data  cards  containing  the 
functional  values  of  the  special  function  must  follow. 

Items  1.2  to  1.30   Format  6X.9I8;  FORTRAN  label  CORBG(I),  1=1,  256 

These  data  cards  are  to  be  supplied  only  if  IDEBUG   ,  ,.  =  0. 
The  values  of  the  special  function  are  to  be  supplied, 
preferably  the  range  of  the  values  should  be  -100,000  to 
100,000,  which  is  internally  transformed  to  the  range  -1  to 
1,  to  make  the  function  comparable  in  norm  to  the  Legendre 
polynomials  and  the  circular  functions. 

F.3.2  Case  Data 

Item  2.1     Format  20A4;  FORTRAN  label  A(I) ,  1=1,  20 

Title  of  the  case.  Columns  77-80  should  be  left  blank  as 
the  text  'BIS'  is  introduced  when  a  portion  of  the  calcu- 
lation is  repeated  after  deviate  points  are  removed. 

Item  2.2     Format  7I3.2D10.0;  FORTRAN  label  N.NPOLY,  NHARJ1,  NPTRE, 
NORRE,  NPREL,  NDUMP,  NDPER,  DFI 
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N   Number  of  data  points.   If  equal  to  zero  computation 

is  stopped. 

NPOLY  Order  of  polynomial  to  be  fitted.   If  IDEBUG   ,„„=  0 
r  '  mod29 

the  highest  order  of  the  polynomial  is  reduced  by  one 

and  a  special  function  is  introduced  instead.   See 

above  for  necessary  change  in  input  data.   Default 

NPOLY  =  3. 
NHAP1!  Number  of  harmonics.   Number  of  associated  linear 

parameters  is  2*NHAKM.  Default  NHARM  =  3. 
NPTRE  Number  of  points  in  a  section  on  which  rejection  of 

erroneous  points  is  enacted.   Default  NPTRE  =  16. 
NORRE  NORRE-1  is  the  order  of  the  approximation  polynomial 

used  in  rejection  of  the  erroneous  points  by  sections. 

Default  NORRE  =  4. 
NPREL  NPREL-1  is  the  order  of  the  polynomial  to  be  subtracted 

from  the  input  data  prior  to  the  estimate  of  the  prominent 

harmonic.   If  NPREL  is  less  than  1  no  action  is  undertaken. 

If  NPREL  is  less  than  zero  the  period  of  the  previous 

case  is  used  as  the  initial  guess.   Subroutines  REJECT 

and  PFRE  are  not  called  in  this  case. 
NDUMP  If  NDUMP  is  greater  than  1  only  points  with  index 

divisible  by  NDUMP  are  retained.  Number  of  points  is 

thus  reduced  to  N/NDUMP. 
DPER  Convergence  criterion  for  the  period,  used  in  subroutine 

XIMIN  (cf.  section  F  2.4).   Default  DPER  =  .05. 
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DFI  Convergence  criterion  for  the  sum  of  squares,  used 

in  subroutine  XIMIN  (cf.  section  F  2.4).   Default  DFI  =  .05. 

Items  2.3  and  subsequent.   Format  6X,  918;  FORTRAN  label  (IY(I),  I  =  1,8), 
NC. 

N  points  of  the  functional  values  at  equidistant  points 
and  the  norm  count.   Norm  count  (of  the  correlator)  is 
subtracted  from  functional  values.   It  should  not  be 
equal  to  zero  as  the  amplitude  of  each  harmonic  is  divided 
by  the  norm  count  for  normalized  amplitude  output . 

F.3.3  Subsequent  Cases 

Items  2.1  through  2.3  are  repeated  for  each  case. 

F.3.4  End  of  Batch 

Items  4.1  and  4.2. 

Two  cards  behind  the  last  case  data  cause  a  normal  stop 

if  the  second  card  has  a  blank  in  the  columns  1  to  3. 

First  card  may  contain  a  text  to  be  printed  before  exit. 

F.4.   FORTRAN  labels  used  in  the  program 

F.4.1  Main  Program 


FORTRAN  Mathematic 

label  symbol  Meaning 

A(20)  case  title  (literal) 

AA  local  variable 
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FORTRAN 
label 


Ma thematic 
symbol 


AMP 

B(15) 
BB(15) 

BIS 

DET 

DFI 

DPER 

F(15) 

FIC 

FIMIN 

FKN 

I 

IC 

IC0RBG(256)  7 

IDBG(20) 

IDEBUG 

IHARM 

IRESP 

IS 

ITIME1 

ITIME2 

IY(256) 

J 

K 

KC 

KK  k 

KS 

N 

NC 

NDUMP 

NGOOD       n 

NHARM 

NORRE 

NPOLY 
NP0LY1 

NPREL 

NPTRE 

NR 

NREJ2 

NREJ4 

PER  t 

PERI  ° 


cb<V 


Meaning 


amplitude  of  prominent  harmonic,  guess  of 

PFRE  subroutine 

parameters  vector 

work  vector  (parameters) 

literal  constant  'BIS' 

dummy  for  MATINV  argument 

convergence  criterion  for  sum  of  squares 

convergence  criterion  for  period 

F-test  values,  Fg5x(i,120) 

confidence  margin  for  sum  of  squares 

minimal  sum  of  squares 

confidence  margin  for  sum  of  squares  factor 

index 

index 

special  function  values 

IDEBUG  options  input 

execution  and  output  options  specifier 

index 

dummy,  response  of  RCHECK  subroutine 

index 

case  execution  start  time 

case  execution  stop  time 

input  data  vector 

index 

number  of  linear  parameters 

index 

number  of  parameters  (including  period) 

index 

number  of  data  points 

norm  count  of  input  data 

data  point  dumping  specification 

number  of  data  points  not  rejected 

number  of  harmonics 

order  1  of  approximation  polynomial  for 

rejection 

order  of  polynomial  to  be  fitted 

number  of  polynomial  columns  in  Z  matrix, 

NPOLY+1 

order-1  of  the  polynomial  to  be  subtracted 

before  calling  PFRE 

number  of  points  per  segment  to  be  examined 

by  REJECT  subroutine 

local  variable 

NPTRE/2 

NPTRE/4 

period,  independent  variable  increment  units 

local  variable 
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FORTRAN 

Mathematic 

label 

symbol 

PH 

PHA 

PHI 

PHIMAX 

PHIMIN 

PHXMA 

PHM1 

RE 

sc 

sec 

SE 

SSR 

TYME 

Y(256) 

YY(3,256) 


YW(256) 

W(50) 

2(256,15) 

ZT 

ZTZ 

ZTZI 


Z,P 

ZT 

Z?ZI 

(z'z)1 


Meaning 

lag,  indep.  var.  Increment  units 

phase  lag,  independent  variable  increment 

units,  guess  by  PFRE  subroutine 

phase  lag,  radians 

phase  lag,  radians,  upper  limit 

phase  lag,  radians,  lower  limit 

lag,  increment  units,  upper  limit 

lag,  increment  units,  lower  limit 

approximate  confidence  region  radius 

deviation  level  for  sum  of  squares  repeat 

criterion  at  residual  examination 

deviation  level  for  point  rejection  at 

residual  examination 

standard  deviation  estimate 

sum  of  squares  for  repeat  criterion 

time  of  case  execution 

data  vector,  good  points  only 

row  1:  data  points;  row  2:  1  for  good 

points,  0  for  rejected  points;  row  3  used 

by  subroutine  PFRE 

work  vector,  residuals 

work  vector 

values  of  component  functions,  matrix  of 

partials 


F.4.2  Subroutine  PFRE 


FORTRAN 
label 


Meaning 


AL  average  minimum 

AH  average  maximum 

AMP  amplitude  estimate 

H(20)  indices  of  maximal  points 

HIST (100)  histogram  of  distances  of  the  extremes 

HMX  local  variable 

I  index 

ID  row  dimension  of  Y 

IH  number  of  maxima  near  peak  in  histogram 

II  index 

IL  number  of  minima  near  peak  in  histogram 

IM  index  of  peak  in  histogram 
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FORTRAN 

label  Meaning 

IMH  local  variable 

IML  local  variable 

L(20)  indices  of  minimal  points 

M3  index  of  row  in  Y  to  be  used  in  subroutine 

N  column  dimension  of  Y 

NH  number  of  maxima  detected 

NL  number  of  minima  detected 

NM  N-l 

NMA  local  variable 

NP  number  of  distances  of  extremes 

PER  period,  independent  variable  increments 

PH  local  variable 

PHA  lag,  independent  variable  increments 

Y(ID,N)  row  1:  functional  values;  row  2:  0  for  rejected  points, 

1  for  retained  points;  row  M3:  used  by  subroutine  for 

calculations . 


F.4.3  Subroutine  ZXIIN  (ZXI) 

Local  variables  and  indices  not  listed. 


FORTRAN 

label  Meaning 

C(256)  cosine  values  vector 

CC  cosine  of  unit  independent  variable  increment 

ISTEP  spacing  of  harmonic  function  evaluations  by  FORTLIB 

function  reference 

N  number  of  original  data  points 

NC  index  of  cosine  column 

NS  index  of  sine  column 

NGOOD  number  of  good  data  points 

NHARM  number  of  harmonics 

NPOLY  number  of  polynomial  columns  -1 

S(256)  sine  values  vector 

SS  sine  of  unit  independent  variable  increment 

XI  period,  independent  variable  increment  units 

XIP  period  at  previous  subroutine  execution 

Z (256,15)  component  function  values  matrix 
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F.4.4  Subroutine  XIMIN 


FORTRAN 

Mathernatic 

label 

symbol 

5(153 

b 

BXC3.15) 

DET 

DFI 

DPER 

ER(N) 

e 

FKH 

FIMAX 

FIMIN 

FINEW 

FIX(I) 

I 

EDEBUG 

II 

IM 

IX 

IW 

J 

K 

(k) 

N 

n 

PER 

t 

R 

0 

X 

XIM 

XINEW 

XINV(3,3) 

W(15) 

WI(15) 

MW(3) 

Y(256) 

y 

Z(256,15) 

z 

Meaning 

parameter  vector 

parameter  vector  storage 

dummy  for  MATINV  argument 

sum  of  squares  convergence  criterion 

period  convergence  criterion 

errors  (residuals)  vector 

sum  of  squares  confidence  margin  factor 

maximal  sum  of  squares 

minimal  sum  of  squares 

new  sum  of  squares 

sum  of  squares  storage 

index 

IDEBUG  options  variable 

index 

maximal  sum  of  squares  point  index 

substituted  point  index 

counter  of  repetitions 

index 

number  of  parameters,  excluding  period 

number  of  rows  in  matrix  Z 

period 

local  variable 

matrix  of  periods,  column  1:  unity  elements; 

column  2:  period  deviations  from  PER  in  ■$ 

P 
units;  column  3:  column  2  squared. 

maximal  period  deviation 

new  value  of  period  deviation  (new  estimate) 

X1 

work  vector,  parameters 

work  vector 

work  vector 

data  points,  good  points  only 

matrix  of  component  functions  values 


F.4.5  Subroutine  LINLSQ 


FORTRAN 
label 


Mathematic 
symbol 


Meaning 


B(K) 

DET 

ERGS) 


parameter  vector 

dummy  for  MATINV  argument 

error  (residuals)  vector 
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FORTRAN 

Ma thematic 

label 

symbol 

FI 

FII 

I 

IDEBUG 

II 

J 

K 

k 

N 

n 

TOO 

Z (256, 15) 

ZPR(15,15) 
ZT(15,256) 
ZTI(15,15) 
ZTZ(15,15) 

y 
z?  i 

(Z^Z)1 
z  z 

F.5.   Output 

Meaning 

sum  of  square  errors 

dummy 

index 

IDEBUG  options  variable 

index 

index 

number  of  parameters  (columns  of  matrix  Z) 

number  of  rows  of  matrix  Z 

data  points  vector 

component  functions  values  matrix 


Output  format  is  illustrated  by  two  sample  problem  outputs.   The 
output  can  be  augmented  by  the  IDEBUG  options  described  in  section 
F.3.1.   These  options  should  be  used  with  discretion  as  a  large  amount  of 
output  can  be  generated. 

F.6.   Sample  Problems 

The  facsimile  of  two  sample  problem  outputs  are  presented  on  pages 
163  to  170. 

F.6.1  Sample  Problem  1. 


The  input  data  is  a  table  of  256  computer  generated  values  of  the 
function  y  =  100,000  cos  (2tt  (i-D/51.3)  +  100,000.   The  values  were 
truncated  to  integers  and  this  truncation  produces  a  certain  noise. 
The  first  parameter,  the  constant  a  in  Eq.  (F-l) ,  is  -.53,  due  to  the 
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truncation,  the  expected  value  being  -.5.   This  value  is  in  the  confi- 
dence interval  of  the  parameter.   The  true  amplitude,  100,000  is  also 
within  the  confidence  interval  of  the  calculated  parameter.   The  slope 
(parameter  2)  and  the  sine  component  are  statistically  insignificant. 

F.6.2  Sample  Problem  2. 

The  analysis  of  the  correlator  output  of  the  test  A15  II,  550  lb, 
correlator  run  7  is  presented.   The  graph  of  the  correlator  output 
(Fig.  34)  displays  two  irregularities  which  can  be  obviously  attributed 
to  the  malfunction  of  the  correlator,  as  the  correlation  function  is 
very  smooth  elsewhere. 

Points  1,  2  and  182  were  rejected  by  the  subroutine  REJECT  prior 
to  the  least  squares  analysis.  Further  points  are  rejected  after 
the  first  least  squares  search,  the  normal  deviate  being  above  5. 
Finally  the  standard  deviation  estimate  is  at  24.12  and  the  residuals 
are  found  to  be  random. 

IDEBUG  was  specified  29  and  the  correlator  response  function  was 
supplied.   It  is  found  to  be  a  significant  component  of  the  correlator 
output  (parameter  5,  value  63.06,  confidence  interval  +  22.56). 
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147 
F  7.   Program  Listing  and  Sample  Outputs 


MAIN 

MAIN  PROGRAM  FOR  ONE  NONLINEAR  PARAMETER  LEAST  SQUARES  ESTIMATE 
IMPLICIT  RE AL*8(A-H,0-Z) 

CI  WENS  I  ON  YY(3,256)  ,  IYI256)  , YI25  6)  ,  YW<256)  ,W(50)  ,Z(  2  56,  15)  , 
1  ZT!15,256),ZTZ(15,15),BB(15),ZTZI(15,15) 
REAL*4  A(20) 
REAL*8  B( 151/15*1.00/ 

F  VALUE  F( I, 120,95 %)     USED,  I  CAN  BE  1  TO  15 
REAL  F(  151/3.92,3.07,2.68,2.45,2.29,2.17,      2.09,2.02,1.96, 

1  1.91,1.87,1.83,1.30,1.77,1.75/ 
COMMON  ZT.ZTZ  ,      ZTZI 
C0MM0N/INT71Y 

CIMENSION  1UBG120) 
COMMON/DBG/  IDEBUG 
READ  9C1, ( IOBGI I ) , 1=1 ,20) 

901  FORMAT (201 3) 
IDEBUG=1 

CO  5  1=1,20 

IF(  IOBGI  I  )  .ECO)  IDEC!I)=1 
5  IDEBUG=IDEBUG*ID8G( I ) 

PRINT  902, I  DEBUG,  ( IDBGt I ) , 1  =  1,20) 

902  FORMAT!'  I DEBUG' , 110,2014 i 

NOTE  NECESSARY  CHANGE  IN  INPUT  IF  MOO ( IDEBUG, 29 ) =0 
DIMENSION  IC0RBG(256) 
IFIMODf  IDEBUG,  29).  ECO)  READ  910,  (  1CCRBG  (  I  )  ,  1  =  1 ,  256  ) 
1  CONTINUE 

CALL  TIME!  I  TIMED 
READ  905, I U I ) , 1=1,20) 

905  FORMAT! 20A41 

PRINT  906, (At  I ), 1  =  1,20) 

906  FORMAT! ' 1' ,20AA,T100, 'PAGE  1') 

READ  9C0,N,NP0LY,NHARM,NPTKE,NCRRr-,NPREL,NDUMP,DPER,DFI 

IF  NPREL  IS  NEGATIVE  THE  PERIOD  OF  THE  PREVIOUS  CASE  IS  USED 

IFtN.EQ.O)  CALL  EXIT 

IF  (NPGLY.CCO)  NP0LY=3 

IF1NHARM.E0.0)  MHARM=3 

IF(NPTRE.ECO)  NPTRE=16 

IF(NORRE.ECO)  NCRRE=4 

IF1DPER.EQ.0.O0)  CPER=.05O0 

IFIDFI.EQ.O. )  DFI=.O5D0 

PRINT  903,N,.NP0LY,NHARM,NPTRC,NCRRE,MPREL,NDUMP,DPER,DFI 

903  FORMAT!'  OPTIONS  SPECIFIED  OR  ASSIGNED  E.Y  DEFAULT,   N=',I4,'    NPO 
1 L  Y  =  ■  ,  I ', ,  •    NHARM=' ,  !'>,'    NPTRE=  '  ,  14  ,  •    N0RRE=',I4,'    NPREL=", 

2  14,'    NDUMP=',I4/'    DPER=" ,F7.3, •    DFI=',F7.3/'  ■) 
900  FORMAT! 713, 2010.0) 

READ  910,  (  !Y(  I  ), 1  =  1, N  1 , NC 

910  F0RNAT(6X,9I8) 

PRINT  911, ( IY( I ),I=1,N) 
PRINT  912, NC 

911  FORMAT!'  INPUT'/'  '/('  ',10110)) 

912  FORMAT!'  NORM  COUNT', 110) 
CO  10  1=1, N 

YY!2,  I  )  =  1. 
10  >Y(1 , I )=DFLOAT! IY( I )-NC) 

IFINPREL.LT. 0)  GOTO  99 
REJECTION  FOR  FIRST   1/4  SEGMENT  AT  95%    CONFIDENCE  LEVEL 

CO  15  I=l,NPTRE 
15    W(I)=YY(1,I) 

CALL  REJECT!  W    , NPTRE ,NORRE , 1 .96D0 ) 
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MAIN 


AT  99%    CONFIDENCE  LEVEL 
.57600) 


CO  16  1=1, NPTRE 
16  YY(2, I )=W( I ) 
NREJ4=NPTREM 
NREJ2=NPTRE/2 

CO  31  I=1,N,NREJ2 
H (I+NPTRE-l.GT.N)  GOTO  31 
CO  25  J=l .NPTRE 
25    HI J)=YY( l.J+I-1) 

REJECTION  FDR  INTERIOR  POINTS 
CALL  REJECT!   U , NPTRE , NORRE t 2 . 
CO  30  J=NREJ4, NPTRE 

30  YY(2, I+J-l )=   W(J) 

31  CONTINUE 

CO  32  1=1, NPTRE 

32  V! I )=YYI l.N-NPTRE+I ) 

CALL  REJECT (W, NPTRE, NORRE, 2. 576U0) 
CO  33  I=NREJ4, NPTRE 

33  YY12.N-NPTRE+I  )  =  M  I  > 
CO  35  1=1,  NPTRE 

35    HI  I  )=YY( l.N-NPTRE+I ) 

RFJFCTICN  FOR  LAST  1/4  SEGMENT  AT  95^  CONFIDENCE  LEVEL 
CALL  REJECT!  W  , NPTRE, NORRE, I . 96DC ) 
NR=NPTRE-NREJ4 
CO  40  I=NR, NPTRE 
40  YYI2.N-NPTRE+1 )=   Will 
PRINT  920 
920  FORMATl'Q  REJECT  SUBROUTINF  CALLEO'/'O') 
K  =  0 

DO  50  I  =  l,.N 
IF(YY(2, I I.EQ.O.DO)  K  =  K  +  1 
50  IFIYYI2, I I.EQ.O.OO)  PRINT  930,1 
930  FORMAT!'  POINT', 14,'    REJECTED') 

IF(K.EC.O)  PRINT  940 
940  FORMAT!'  NIC  POINTS  REJECTED'  ) 

PFRE  CALLED  MOW  WITH  OPTION  OF  FIRST  SUBTRACTING  A  POLYNOMIAL 
IFIMPREL.LT. 0)  GOTO  99 
IFINPREL.LT.2)  GCIO  90 
CO  60  1=1, N 
60  Y( I )=YY( 1,1) 

STORE  FUNCTIONAL  VALUES 
DO  70  1=1, N 
Z(1i  1>  =  1.D0 
CO  70  J=2,NPREL 
70  7! I, J)=Z! I , J-1)*I 

CO  75  1=1,15 
75  81! 1=1.00 

CALL  WMTRA!Z,ZT,N,NPREL,256,15) 
CALL  V.;mPRD(ZI,7,  7TZ.NPREL.N,  15,256) 
CALL  MATIMV(ZTZ,MPREL,B,1,DET,15) 
CALL  WWPRDV(ZT,Y,W, NPREL.M, 15  1 
CALL  V.'MPROV!  ZTZ,'.'. ,  B,  NPREL  ,  NPREL  ,  15  ) 
CALL  WMPROVIZ  ,e,YV.,N,NPREL,256) 
CO  80  1  =  1, N 
80  YY(  1,1  )=Y!  D-YWI  I  ) 
90  CALL  PFRE(YY,3,N,3,AMP,PER,PHA) 
PRINT  950, AMP, PER, PHA 
950  FORMAT! 'OINI HAL  GUESS ' / • OAMPL 1 1 UCE • , T 1 8 , FB . 2/ '  PER  100' , T 1 8 , Fl 0. 4/ 
1  '  LAG' , • (PHASE) ' ,T18,F10.4) 
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RETAINED' ) 


THE  NORMAL 

■  0 


MAIN 

IFINPREL.LT. 2)  GOTO  100 
CO  95  1=1, N 
95  YYU,  I  )  =  Y(  I  ) 

GOTO  100 
99  PRINT  955, PER 
950  FORMAT  I 'OPERIOO  OF  PREVIOUS  CASE  USED  FOR  INITIAL  GUESS,    PERIOD" 

1,1  1  0.  4  ) 
100  K=l 
C      DUMP  POINTS  ACCORDING  TO  SPECIFIED  NDUMP 
IF  (NDUMP.LE.O)  GOTO  106 
CO  105  1=1 ,N 

105  IF1M0DI I,NDUHP).NE.O)  YY(2,I)=0.D0 
PRINT  956,N0UMP 

956  FORMATt'OONLY  POINTS  WITH  INDEX  DIVISIBLE  BY1   ,13," 

106  CONTINUE 
C       INITIALIZE  VECTOR  Y  AND  ARRAY  Z  WITH  GOOD  POINTS  ONLY 

CO  110  1=1, N 

Y(K)=YY(1, I ) 

Z(K,2)=DFL0AT( I -N/2) /IN/2) 
C      FOLLOWING  STATEMENT  INITIALIZES  COLUMN  NP0LY1  OF   Z   TO 
C      VALUES  OF  THE  CORRELATCR  BACKGROUND,  IF   MOD ( I  DEBUG , 29 ) 
C      NOTE  NECESSARY  CHANGE  IN  INPUT 

110  lFi?Yi?I:??^:g?^CK2K  +  lZ,K,'P0LY+11=  IC0RBG""^ 
NG00D=K-1 

IFINGOOD.EQ.O)  GOTO  172 
NP0LY1=NP0LY+1 

IF ( MOD U DEBUG, 29). EQ. 01  NP0LY1=NP0LY 
CO  120  I=1,NG00D 
Z(  I, 11  =  1. DO 

ifinpolyi.lt. 3)  core  119 

DO  110  J=3,NPCLY1 

lie  ni,j)=( i2*j-i)*z(i,j-i)*z(i,2)-(j-i)*zii,j-2))/j 

120  CONTINUE 

IF(«00( ICEBUG.2 l.NE.O)  GOTO  124 
DO  123  I=1,NGCCD 

123  PRINT  958, I, (Z( I , J ) , J= I , NP0LY1 ) 

124  CONTINUE 

95B  FORMAT!'  • , I  3, ( T5, 10F10.5) ) 
C      INITIALIZE  SUBROUTINES   ZXI 
CALL  ZXI  IN(Z,NPOLY,NHARM,N) 
K*NP0LY+2*NHARM+1 

FKN=K*F(KK)/(NGCOD-KK) 

CALL  XIMIN(Z,a,Y,YW,NGCOD,KfPER,DPER,FIMIN,DFI,FKN) 

PRINT  960!NPOLYlNHARM,(I,B(l) , I = 1 , K ) ,KK, PER 
FORMAT  (  •  1'  ,20 A'.,  T  100,  'PAGE  2'  I 


959 

960 


FOU'ATI  '(MINIMUM  FOUND •/• OORDER  CF 
1ARH0NICS',  I  3  /^PARAMETERS'/'  '/I' 
CO  121  1  =  1,  N 

121  Y(I)=O.DO 

CO  122  I=1,NGC0D 

122  YlIDINTIZ(I,2)*N/2+.l  +N/2))  =YK(I) 
PRINT  961, (Y( I),I=1,N) 

961  FORMATCORESICUALS   (EXACT  ZERO  CFNCFES 
1C  '.1CF10.2)) 


POLYNOMIAL'  ,13, 
' ,I3,G15.6) ) 


NUMBER  OF  H 


REJECTED  POINTS)1/'  •/ 
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MAIN 

IF(MOD(IDERUG,23).NE.O)  GOTO  117 
PUNCH  975, (Y( I) , 1=1, N) 
975  FCRMAT(6X,9F8.2) 
117  CONTINUE 

SE  =  DSQU(FIMIN/(NCCOD-K-l)  ) 
FIf>FKN*FIMIN 

PRINT  966,SE,F!MIN,FIC  „„  . 

966  FORMAT!  'OSTANDARO  DEVIATION  ESTIMATE'  ,F8.2»«         SUM  OF  SQUARES 
1=',F16.0,'    95";  CONFIDENCE  MARG  I  N=  ■  ,  F  16  .  0  ) 
CALL  RCHECK(YW,NGOOD,IRESPl 
1?6  CONTINUE 
C      DETERMINING  THE  LINEARIZED  CONFIDENCE  LIMITS 

C      IF  NUMBER  OF  PARAMETERS  IS  GREAIER  THAN  FIFTEEN  ONLY  ESTIMATE  OF 
C       LINFAR  PARAMETERS  IS  AVAILABLE 

PRINT  964, (At  I ) ,1  =  1,20)  .  ,   . , 

964  FORMAT! '1' ,20A4,T100, 'PAGE  3 • / • OL I  NEAR  I  ZED  CONFIDENCE  LIMITS'/'  ') 
IFINP0LY*2*NHARM+2  .GT.  15)  PRINT  965 

IF(NPCLY+2*NHARM+2  .GT.  15)  KK  =  K 
IF(NP0LY+2*NHARM+2  .GT.  15)  GOTO  135 

965  FORMAT!'  TCO  MANY  PARAMETERS-PERIOD  CAN  NOT  BE  INCLUDED  IN  THE  CON 
1FIDENCE  LIMITS  CALCULATION') 

B(KK)=PER 

PER  1 =-2. 00*3.141 5926/ PER**2 
DO  130  I=1,N3C00 
71 I,KK)=O.UO 
CO  130  J=1,NHARM 
KC=NP0LY+2*J 
KS=NP0LY+2*J+1 
130  /( [,KK)=Z( I ,KK)+PER1*(-B(KC)*Z(I,KS)+B(KS)*Z(I,KC) )*(Z< I,2)*N/2+ 

1  N/2-1) 
IFIMODt IDEBUG.31.NE.0)  GOTO  134 
CO  133  I=1,NGC0D 

133  PRINT  958, I, (Z( I,J),J=1,KK) 

134  CONTINUE 

CALL  WMTRAU.  ZT  ,  NGOOD  ,  KK  ,  256  ,  15  ) 

CALL  WMPR01ZI , Z,7TZ,KX,N GOOD, KK, 15,256) 

CALL  MATINVIZrZ.KK.W, 1.DET.15) 

135  CONTINUE 

CO  138  1=1, KK 
133  'MI  )  =  DSCRT  (DABS(FKN*FIMIN*ZTZ(  1,1))) 

PRINT  9  7C, ( I,E! I ) . Wl I ) , I=1,KK) 
970  FORMAT!"  ',I3,F15.2,'     +/-',G15.5) 

P,CAL*4  BIS/' BIS  •/ 

SSR=0.D0 

SC=3.291*SE 

SCC=3.09*SE 

CO  125  1=1, N 

IFIOAPSI  Y(  I  )  )  .GT.  SO  SSR  =  SSi!  +  Y  (  I  )  **2 

IFIOABSIY! I ) )  .GT.  SCO  YY(2,I)=C.D0 
125  IFIOABSIY! I ) )  . GT .  SCO  PRINT  962,1 
962  FORMAT! '  POINT', 14,*   SHOULC  BE  REJECTED') 

IF1SSR.  LT.  FIC/2)  GOTO  139 

PRINT  963 
961  FORMAT!'  CALCULATIONS  SHALL  BE  REPEATED  WITH  MORE  POINTS  REJECTED' 
1  ) 

A(19)=BIS 

K  =  l 

ccro  IC6 


151 


MAIN 

139  CONTINUE 

DO  149  IHARM=1,NHARM 
IC=NP0LY+2*IHARM 
tS=NP0LY+2*IHARM+l 
AMPL  =  DSQRr(B( IC )**2+8 1 I 5 )**2) 
RE=DMAXl(DABSJWIIC)  ),DABS(W(  IS))) 
^'Rf;GT.AKPL)  PRINT  805,  HARK 
IFIRF.GT.AMPL)  CfirO  148 

8°'J  pTK^BIIcI^fl?,.™*  H«"«W*C.I3..    NOT  S.GNIF.CANT., 
OPHI WANGLE (AMPL, RE) 
PHIMAX=PHl+DPHI 
PHIMIN*PHI-DPHI 
ANPN=AKPL/NC 
»A=3.1415926*2/PER 
PH=PHI/AA 
PHMA=PHIMAX/AA 
PHMI*PHIMIN/AA 


149  CONTINUE 

DO  150  1=1, KK 

150  iUI!j8s?Rr(PAtS(i'TZ"'I>>> 

CO  160  1=1, KK 
CO  160  J=1,KK 

PRlf'T,320^rZ1I*J,/IW,I,*H(J,, 

820  D0R170,I°1AKKMETERS  CnRRELATIOrJ  f'ATRIX'/'O') 

8372CFC^LAT{;NEi!?!^!°G1?-5/'  '."0.5GX2.5I 

TYME=( ITIME1-ITIME2)*I-.1D-1) 

PRINT  840.TYME 
840  FORHATCOCASE  COMPUTATION  TIKE  i.F10.2,«   SECONDS') 

END 


SUDR 
IMPL 
DIME 
CO  MM 
CALL 
COHM 
I F  ( M 
CO  1 
1  PRIM 

11  CONT 
120  FORM 

CO  4 
DO  4 
ZTZ  ( 
CO  4 

A  ZTZI 
DC  6 
CO  6 

6  Z T  Z  ( 
DIME 
CALL 
I F  ( M 
DO  3 

3  PR  IN 

12  CONT 
CALL 
CALL 
I  F  (  M 
CALL 
DO  5 

5  PRIN 

13  CONT 
100  FORM 

CALL 
CO  1 
10  FR(I 
I F  (  M 
PRIN 

19  CONT 
FI  =  0 
CO  2 

20  FI  =  F 
F  I  I  = 

21  CONT 
RETU 
FND 


OUT 
IC1 
NSI 

ON 

V.f 
ON/ 
0D( 

1  = 
T  1 
INU 
ATI 

1  = 

J  = 
I.J 

1  I 
I.J 

1  = 

J  = 
I.J 
NSI 

WM 

om 
1  = 

T  I 

INU 

KM 

MA 

0D( 

WM 

1  = 

T  1 

INU 

ATI 

WM 

0  T 

)  =  Y 

out 

T  1 

INU 
.DO 
0  I 

i+rj 

FI 

INU 
RN 


INE  LI 

T  REAL 

ON 

ZT.ZTZ 

TRAIZ, 

DBG/  I 

IDEBUG 

l.N 

20,1,1 

E 

.13 
l.K 
I.K 

1=0. DC 
=  1,N 
)=ZTZ( 
l.K 
1,1 

)=ZTZ( 
ON  ZTI 
ASGNIZ 
IDEBUG 
l.K 
20,1, ( 

PRDVIZ 

T  I  NV  I  Z 

IDEBUG 

PRUIZT 

l.K 

00, (ZP 

E 

•  ',10 

PR0VIZ 

=  1,N 

I  I  l-ER 

IDEBUG 

00, IER 

E 

=  l,N 

ad  )**; 


LINLSQ 

NLSQ(Z,B,Y,ER,N',K,FI) 
*8(A-H,0-Z) 

B(K), YIN), ER(N),ZT<15, 256), ZTZU5, 15), 7(256, 


ZT,N,K,256, 15) 

DEBUG 

,5). NE.O)  GGTC  11 

Z(  1,J),J=1,K) 

, (10G12.3I ) 


I,J)+ZT( 1,1 l)*Z(II,J) 


152 

56,15) 

J.  I) 

115, 
TI  ,Z 
,7). 


15) ,ZPRI15,15) 
TZtK.Kt 15) 

NE.O)  GOTO  12 


ZTI  (  I , J) , J  =  1,K) 

T,Y,B,K,N,15) 
TZ,K,B,1,DET,15) 
,11)  .NE.O)  GOTO  13 
I,ZTZ,ZPR,K,K,K, 15, 15) 

K( I, J) ,J=1,K) 


G12. 
,B,E 

(  I  ) 

,13) 

(I), 


A) 
,N,K,256) 


.NE.O)  GOTO  19 
I=1,N) 


153 


100 


10 
15 


ZXIIN 

SUgRClUTI*JK    ZXIINU.NPOLY.NHAUf'.N) 
IMPLICIT ■ REAL*8(A-H,0-Z) 
XIpIo^O0N    Z("6'l'J,'S(256),C(256) 
If  ^?C)Vy  +  2*r'JHARr-,+  1-GT.  15)     PRINT     100 

i?MMi;*3SRBISHSfJf?}grN5§  T0D  LARGE*  QUIT  CASE,) 

PI  =  3. 141592653 

RETURN 

ENTRY  ZXKZ.XI.NGOQD) 

IF(XI.EQ.XIP)  GOTO  15 

XIP=X1 

ISTEP=10 

DO  10  I=1,NHARH 

ARG=2*PI*l/XI 

CC-DCOS(ARG) 

SS=DSINIARG1 

ISYEP=10 

CO  5  n*l,Ni  ISTEP 

IE=I I+ISTEP-1 

IFUE.GT.N)  IE=N 

CIII  l«DCOS(ARG*l II-l) ) 

S(II)=OSIN{ARG*(n-l)) 

112=11+1 

DO  5  J=fl2.IE 

SU)  =  S(J-1)*CC  +  CU-1)*SS 

C<J)=C(J-1)*CC-S(J-1)*SS 

NS=NP0LY+l+2*J 

NC=NPCLY+I*2 

CO  10  J=1,NG00D 

INDX=IDINT(Z(J.2)i 

ZUtNC)=C< INOX) 

Z(J,NS)=S( INDX) 

RETURN 

END 


N/2+N/2+.1D0) 
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XIMIN 

IDi5Hw{l5)fwi?ll)5^!i?,,Y(N)'ER{N,,XINV.(3»3,'X(3'3J'FIx«3>tBXt3, 

SRJIK8IMI13J?    ?INV,X,FIX,8X,H,WI  ,HH 

COMMOM/DBG/  IOEBUG 

CO  10  1=1,3 

X!  I,1)=1.D0 

XI  I, 2)=-  2. 00  +  1*1. DO 

X(  1,3)=X(  1,2)**2 

CALL  ZX1  (Z.PER  +  X!  I,2)*CPER,\|) 

CALL  LINLSCIZ.W, Y.ER.N.K, FIX! [)) 

LU  1  U  J-lfK 

10  BX(I,J)=U(J) 

FlMIN  =  Di,INl(FIX!l),FIX(2),FIX(3)  I 
12  R=  FIXI3)-FIX(2))/(FIX(2)-F1X(1 

'F'R;LT..95D0  .OR.  R.ST. 1.0500)  GOTO  50 

UU  Yj     \  —  I  (  K 
14  W( [)=BX(2, 1 ) 

rx=i 

IF(FIX(3).LT.FIX(1) )  IX=3 

FIX(2)=FIX( IX) 

X(2,2)=X( IX, 2) 

X(2,3)=X( IX, 3) 

CO  17  1=1,15 

BX(2, I )  =  3X( IX, I  ) 

IF1IX.F0.3)  X(3,2)=  2*X(2,?)-X( 1,2) 

CALL  ZXI  (Z,PEK  +  X( IX,2)*DPES,N) 

CO  18  j=lLKC^'rt'Y*ER'N'K*FIX(rXn 

EX( IX, J)=w( J) 

GOTO  12 

CALL  WMASGN(XINV,X,3,3,3) 

V.  W  (  I  )  =  F  I X  (  I  ) 

CALL  MATINV{XINV,3,WH,1,DET.3] 

XINEW=-.5*WH(2)/WW(3) 

IW  =  0 

CALL  ZXI (Z,PER+XINEW*DPER,N) 

CALL  LINLSQ(Z,V;,Y,ER,\,K,FINEW) 

ln^"r     Pv^0,'1?1?!0-01    Pimr    101.XINFW.FINEW 

FORMAT!'  XINbW=«,F8.4, '    F INEW= • , G 1 5 . 6 ) 

FIKIN=0MIH1  FIMlft.FINEW) 

FIMAX=FIMIN 

IH*0 

CO  60  1=1,3 

FIMAX=0MAX1(FIMAX,FIX(I)J 
FfFIXm.EQ.FIHAX)  [l'=I 
IFUM.EQ.O)  PRINT  100 

XIMIN,  AT  60') 


17 


18 

50 


51 


101 


60 

100 


FORMAT (•  ERRCR 
IFUM.EQ.O)  IN 
CO  65  1=1, K 
65  hi! I )=8X(IM,I I 
FIX! IM)=FINEw 
X! IM,2)=XINEK 
XUM,3)*XINEW**2 


IN 
1 
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66 
70 


135 


74 

76 


110 


72 
73 


71 


80 


X  I  M 
CO  6 
[F(F 
CD  7 
P.X(  I 
I  F  (  M 
IF(D 

1  1.0 
IF(F 
1 1=1 
[F(  [ 
IF(  I 
FORM 

IHICH 
IF(X 
IF(X 

co' il- 
ea 7 

X(  [  , 

PRIM 

FORM 

1C    P 

goto 

XINE 
GOTO 

XI  fJE 
GOTO 
CGNT 
IL=l 
DO    8 

i  r  (  f 

PER  = 

00 
B(l) 
RETU 
ENO 


XIMIN 

=  XU,2)  ■  ' 

6  1=1,3 

tX(I).EQ.FlMIN)   X1M   =  XU,2) 
0  1  =  1, K 

y,  n=wi  i ) 

00  {  10EBUG.19)  .ECO)  GCFC  T>  ,,  „. 

AX1<X<1,2),X(2,2),X<3,2)>-DVIN1<X(1,2),X(2,2),X(3,2)).GT. 

0)  GOTO  50 

IMAX/FIMIN  .GT.DFI*FK'i  +  l.C0)   GCTC  50 

W.GT.3)  PRINT  135 

AT('""xiMIN  FORCED  EXIT,  XH'IN  IS  NOT  BETWEEN  TWO  POINTS  U 

ER  SUM  OF  SQUARES' ) 

IK.EQ.   0FIN1(X(1,2),X(2,2),X(3,2) ) 1G0T0  72 

IM.EQ.   0rAXl(X(l,2),X(2,2),X(3,2)))GOTO  73 

DUE 

6  1=1,3 

2)=PER+X( I ,2)*DPER 

T  110, (X(I,2) ,FIX( I ) ,  1  =  1,3) 

ATI'O  SUBROUTINE  XIMIN  EXIT     LAST  THREE  POINTS  ..'/ 

ERIOD='  ,F8.4,'    SUM  OF  SCUARES=  •  ,  G  1  5.  5 ) ) 

71 
W  =  Xr<-1.2-0t-'AXl(X(l,2),X(2,2),X(3,2)  )  *.2 

55 
W=XIM*1.2-DHIN1 (X(1,2),X(2,2),X(3,2))  *.2 

55 
INUE 

0  1=1,3 

IXdl.EQ.FIHIN)  IL=I 

X  (  I  L  ,  2  ) 

£5  1=1, K 

=BX(IL,I) 

RN 


TH 
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10 


11 


13 

15 


18 


20 


25 


30 


SUB 
INI 

IMP 

DIM 

CON 

PC 

PHA 

AMP 

M',= 

CO 

L(l 

HI  I 

SHI 

FRR 

CO 

YIM 

IF( 

NMA 

NL  = 

M!  = 

FIN 

CO 

IF! 

IYIM 
IF( 

1Y|M 
GOT 
NH= 
HI 
GOI 
M_  = 
LIN 
CO 'J 
CON 
DO 
HIS 
CO 

II 

HIS 

CO 

II 

HIS 

FMX 

If  = 

CO 
IFI 

HMX 

per 
np 

[ML 

fMH 

CO 
PER 
NP  = 
IFI 
IFI 


PFRE 

ROUTINE  PFRE(Y,ID,N  , M3 , AMP , PER , PHA ) 
TIAL  ESTIMATE  OF  THE  PROMINENT  FREQUENCY 


LICIT 
ENSIO 
MON/I 
=  0. 
=  0. 
=  0. 
N-l 

1  1  =  1 
|=0 

)  =  0 
FT  FU 
ATIC 
10  1  = 
3.1)  = 
YI2,  I 
NM-4 
0 
0 

D  THE 
15  1  = 

Y  ( M  3 , 
3,  1+1 

Y  I  M  3  , 
3,  1  +  1 
0  15 
NH+1 
H)=I 
0  15 
NL  +  1 
L)  =  ! 
TINUE 
STRUC 
18  1  = 
T(  I  )  = 
20  1  = 
H  1 1  )  - 
Till) 

2  5  1  = 
LI  I  )- 
Till) 
=  HIST 
1 


REAL*8IA-G,C-Z) , INTEGER  (H) 
■1    YUD,  N)  ,L(?0)  ,H(20)  .HISTUOO)  ,ICUMY(116) 
MT/  HIST.L.H, IDUMY 


,20 

MOTION  FROM  ROW  1  TO  ROW  M*3  AND  CONSTRUCT  SUBSTITUTION  FOR 

POINTS 

2,  MM 

Yd,  I  ) 

). EC 0.)Y(M3,I)=(Y(lf  I- i)  +  YU, 1  +  1)1/2 


LOCAL  MAXIMAL  AND  MINIMAL  POINTS 
6,  NMA 

I) .GT.0MAX1 (YIM3, 1-4 ) , Y I M3 , 1-3 ) ,Y(M3,I-2) ,Y(M3, 1-1) , 
) ,Y(M3, 1+2) ,Y(M3, I+3I.YIM3, 1+4) I )  GOTO  11 
I I.LT.DMIN1 !YIM3,l-4),Y(M3, 1-3) ,Y(M3, I-21.YIM3, 1-1 ) , 
),Y(M3,I+2)fY(M3,l+3)tY(M3fI+4) >)  GOTO  13 


T  HISTOGRAM  OF  THE  DISTANCES  BETWEEN  ADJACENT  MAXIMA 

I,  100 

0. 

2,'JH 

11(1-1) 

=HISII  ID  +  1 

2,NL 

LI  1-1) 

=  HIS)  I  II1+1 

I  1) 


1=1, 100 
X.LE.HIST1  I  )  )  IM: 
AXO(HMX,HISTl I ) ) 


2 
M  +  2 

1  =  1  ML, IMH 
ER+I*HIST( I ) 
+  HIST I  I ) 
.LT.l )  GOTO  60 
.LT.l)  PER=0. 
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PFRE 


40 


50 


55 


60 


PER 

CO 
IF( 
DO 
IF( 

in 

IF( 
IH= 

IL  = 

PH  = 
A 11= 
AL  = 

Co 

IF( 
PH= 
IH= 

AH  = 
CON 
CO 
IF( 
PH  = 
IL  = 
AL  = 

cc;»j 

PHA 
FOL 
PHA 
AMP 

RET 

END 


PER 
'.0  I 
DABS 
45  I 
DABS 
DABS 
DABS 
0 
0 
0. 
0. 
0. 

50  I 
HI  I) 
PH  +  D 
IH  +  1 
All  +  Y 
TINU 
5  5  I 
LI  I) 
PH+D 
IL+1 
AL  +  Y 
TINU 

=  PH 
LOW! 
=  6.2 

I  AH 
URN 


/NP  ■ 

=  2,MH 

(HI  I )-H(I-ll-PER).CT.2.D0)  H(I-1)  =  0 

=  2,NL 

(L( I )-L( I-1)-PER).GT.2.D0)  LI  1-11=0 

(H(NH)-H(NH-l  1-PER)  .GT.2.C0)  H(NH)  =  0 

(LI  NL)-L(NL-1)-PER).GT.2.D0)  L(NL)=0 


=  1,NH 

.EO.O)  GOTO  50 

(■lODIDFLOATIHI  I  )  l.PER) 

1 1-'  3  ,  H  I  I  )  ) 
E 

l.NL 

EQ.O)  GOTO  55 
HODIDFLCATILI I ) ) +PER/2 , PER ) 

(M3,L(U) 

E 

/( IH+IL)  -I. 

NG  STATEMENT  EXPRESSES  PHASE  LAG  IN  RADIANS  WHEN  C  IS  REMOV 

83185*PHA  /PER 

/IH-AL/ID/2. 


ED 
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RCHECK 

SUBROUTINE  RCHECK ( Y,M, IRE SP  ) 
CHECK  PERFORMED  ON  THE  NUMBER  OF  RUNS  OF  THE  SAME  SIGN 

IMPLICIT  REAL  *8(A-H,0-ZJ 

DIMENSION  YIN)iK(256I 

COMMON  /INT/  K 

NRUN=1 

NMM=N-T 

CO  10  1=2, N 
10  IF(Y( I )*Y( 1-1 I.LT.O.)  NRUN=NRUN+1 

NP  =  0 

NM=0 

CO  20  1=1, N 

IF(Y(   I  I.LE.O.  )  NM=NM+1 
20  IF(Y(   Il.GT.O.)  NP=NP+l 

MI=2*NP*NM/ INP+NM )+l 

SIG=DSQRT(2.C0* 

1  NP*NM*t2.*NP*NW-NP-NM)/t ( NP  +  NM 1**2* (NP  +  NM- 1. )) ) 
Z=!NRUN-MI+.5DC)/S!G 

IF(DABS!Z!.LT.1.645DC)  IRESP=0 
IF(DARSIZ)  .LT.1.64r>D0)  GOTO  30 
IFIZ.GT.O.)  IRESP=1 
iriZ.LT.O.I  IRCSP=-1 
PRINT  100,MI,NRUN,Z 
CO  2  5  1=1 ,N 
K( I  1=0 
?5  IF(Y(   I I.GT.O. )K( I )=1 
PRINT  110, (K( I ) , 1=1 ,NI 
100  FORMAT! 'OSEQUENCF  OF  SIGNS  IS  NOT  RANDOM  AT  THE  90^  SIGNIFICANCE  L 
1FVEL,  EXPECTED  NO.  OF  RUNS  IS', 14/'    ACTUAL  NO.  I  S '  ,  1 4 ,  '    NORMAL 

2  DEVIATE  IS',F8.3  /'OSEOUENCE  OF  SIGNS'/'O1) 
110  FORMAT! •  ',12811) 

GOTO  AC 
30  PRINT  120,MI,NRUN,Z 
120  FORMAT! 'GRANDCMNESS  NOT  REJECTEO  AT  ^0":  CONFIDENCE  LEVEL,  EXPECTED 
1  NO. OF  RUNS  IS", 14,'    ACTUAL  NO.  OF  RUNS  IS', 14,'    NORMAL  DEVIAT 
2E  IS' ,F8.3) 
40  RETURN 
END 
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REJECT 

SUBROUTINE  REJECT! Y, H,N, ESS ) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  Y<M),YY(5G),X(50,1C),XT(10,50),XTX(10,10),XTXI<10,1 
lREAL*0  B(10,/10,IJ!^J01•WWl50*5C,  ^"UHYC36I 
COMMON  YY,X,XTX,WW 
COMMON/INT/XTXI,  IOUf-'Y 
COMMON  /A/  W.XT 
CATA  HP/O/tNP/O/ 

IF(M.EC.MP.AND.N.EQ.NP)  GOTO  40 
CO  10  1=1, M 
XII, 1)=1.D0 
CO  10  J=2,N 
10  X( I, J)=X< I , J-l)*[ 
CALL  WMTRA(X,XT,M, 


0), 


20 


40 


44 
45 


50 


.NtSO.10) 
WMPP,U(Xr,X,XTX,N,M,N,  10,50) 
WMASGN(Xrxl,XTX,N,N.lO) 
MATINV(XTXI,N,B,l,0ET,10) 
WMPRD(XTXI,XT,W,N,N,M, 10,  10) 
WMPRD  I X  t h  «UV< ,  K, N,  Mf  50,  10 ) 
I  - 1 ,  M 


CALL 

CALL 

CALL 

CALL 

CALL 

CO  20 

KW (  I,  I  )  =  WW(  I,  I  1-1. 

CALL  WMPRD (XT  XI ,XTX,W,N,N,N, 10,10) 

fvP  =  N 

MP  =  M 

SS  =  0. 

CO  45  1=1, M 

YY( I )  =  C. 

CO  44  K=1,M 

YY( I )=YY! I  )+,-:w(  I,K)*Y(K) 

SS=SS+YY( [ )**2 

S=DSQRT(SS/(M-N) ) 

CO  50  1  =  1, M 

YII 1=1. DO 

IHDABS(YY(I)).GT.FSS     *S)  Y(I)=O.DO 

END 
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WHAT 

SUBROUTINE  InMAT  (  A  ,  B  ,  C  ,  AA  ,AT  ,  N,N,  I ,  NFUL  ,  MFUL  ,  V  ,  VI ) 

REAL*  8  A (NFUL, «  ),BIHFUL,L  )  ,  C  (  NFUL  ,  L  )  ,  A  A  I  NFUL  ,.M  >  ,  AT  <  KFUL  ,  N  )  ,  V  (  H  ) . 
1  Vi(N) 

ENTRY  KMPRDIAf B,CtN,M,L,NFUL,MFUL) 

CFJ  10  1=  1 ,  N 

CO  10  J=1,L 

C( I, J)=0. 

CO  10  KK=l,M 
10  CI  I, J)=CU, J )+A( I,KK)*B(KK, J) 

RETURN 

ENTRY  Wiv,TRA(A,AT,N,M,NFUL,N'FUL) 

CO  20  1  =  1,  N 

CO  20  J=1,M 
20  ATI  J,  [  )  =  A( I , J) 

RETURN 

ENTRY  WMASGNIA,AA,N,K,NFUL) 

CO  30  1=1,  N 

DO  30  J=1,N 
30  H I, J)=AA< I , J) 

RETURN 

ENTRY  KM ADD ( A , AA, N, M, NFUL ) 

CO  40  1=1, K 

DO  4  0  J=1,M 
40  A( 1, J)  =  A( I, JI+AAI  I, J  ) 

RETURN 

ENTRY  HMSUB(A,AA,N,M,NFUL) 

CO  50  1=1, N 

CO  50  J  =  l,f" 
5U  A( I, J)=AA< I, J)-A( I, J) 
RETURN 

ENTRY  WMNUl  (  A  ,N  ,  f' ,NFUL  ) 

CO  60  1  =  1,-1 

CO  60  J=1,M 
60  All, J)=0. 

RETURN 

ENTRY  WMUN I T I  A , M , M ,NFUL ) 

K=HINO(NFUL,M) 

CO  69  1=1, N 

CO  69  J=1,M 

69  A(I,J)=0. 
CO  70  1=1, K 

70  All,  I  1  =  1. DO 
RETURN 

ENTRY  WMPR0V(A,V,K,N,M,NFUL1 
CO  80  1=1, N 

1  mu«o, 

CO  RO  K=1,M 
80  MI  )=H(I  )  +  A(  I  ,K)*V[K) 

RETURN 

ENTRY  WVPRCM(W,A,V,N»M,NFUL) 

DO  90  1  =  1,  C. 

V( I )=0. 

CO  90  K=1,N 
90  VI  I )=V( I l+W (K)*A(K,  I ) 

RFTURN 

END 


161 


HAT INV 

SUBROUTINE  MATINV  I A,N,B,M,DETERf> ,NMAX.) 

C      KATRIX  INVERSION  WITH  ACCOMPANYING  SOLUTION  OF  LINEAR  EQUATIONS 

C****  REMOVE  NEXT  STATEMENT  IN  SINGLE  PRECISION  VERSICN 
IMPLICIT  REAL*8  (A-H.O-Z) 
REAL*'.  PIVOT 

DIMENSION  A(NMAXtN).  BINMAX.N) 
COMMON  /F402/  PIVOT (100),  INOEX(IOO) 

C      INITIALIZE  DETERMINANT  AND  PIVOT  ELEMENT  ARRAY 


C 


DETERM=1.0 
CO  20  1=1, N 
PIVOT! I  1  =  0.0 
20  CONTINUE 
C 

PERFORM  SUCCESSIVE  PIVOT  OPERATIONS  (GRAND  LOOP) 

CO  550  1=1, N 
C 

SEARCH  FOR  PIVOT  ELEMENT  AND  EXTEND  DETERMINANT  PARTIAL  PRODUCT 

AMAX=0.0 

CO  105  J=1,N 

IF  (PIVOT! J). ME. 0.0)  GO  TO  1C5 

CO  100  K=1,N 

IF  (PIVOT(K).ME.O.O)  GO  TO  100 

TENP=DABSI A( J,K) ) 

IF  (TEMP.LT.AMAX)  GO  TO  100 

IROW-J 

ICOLUM=K 

AMAX=TFMP 
100  CONTINUE 
105  CONTINUE 

INDEX( I l=4096*IRCW+ICCLUH 

J=IROW 

*MAX=A( J, ICOLUM) 

CETERM=AMAX*DETERM 

RETURN  IF  MATKIX  IS  SINGULAR  (ZERO  PIVOT)  AFTER  COLUMN  INTERCHANGE 
IF  (DETERM.EC.O.O)  GO  TO  600 


C 


c 


PIVOT! ICOLUM)=AMAX 


C      INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

IF  ( IRCU.EQ. ICOLUM)  GO  TC  260 
CETERM=-DETERM 
DO  2  00  K=l,N 
SKAP=At J,K) 
A(J,K)=A1 ICOLUM, K) 
A( ICOLUM, K)=SWAP 
200  CONTINUE 

IF  (M.LE.O)  GO  TC  260 
CO  2  50  K=1,M 
SWAP=BIJ,K) 
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MATINV 

B(,),K)  =  B(  ICCLUM.K) 
B( ICOLUM,K)=$WAP 
250  CONTINUE 
C 

C      DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 
C 

260  K=ICOLUM 

A(  TCPLUM,K)  =  l.O' 
DO  350  K=1,N 

A(ICPLUM,K)=A(ICOLUM,K)/AMAX 
350  CONTINUE 

IF  (M.LE.O)  GO  TO  380 
00  370  K=1,M 

B(  ICOLUM,K)=B(  IC0LUM.K1/AMAX 
370  CONT IMUE 
C 

C      REDUCE  NCN-PIVCT  ROUS 
C 

380  CO  550  J=l,N 

IF  (J.EQ.ICOLUK)  GO  TO  550 
T  =  Al  J.ICOLUM) 
A(  J,  I  COLO.")  =0.0 
00  450  K=1,N 

A(  J,K)=A(  J,K)-A(ICOLUM,K)*T 
450  CONTINUE 

IF  (M.LE.O)  GO  TO  550 
BC  500  K=1,M 

P(  J,K)=B(  J.KJ-BI ICOLUM,K)*T 
500  CONT  I'lUE 
550  CONTINUE 
C 

C      INTERCHANGE  COLUMNS  AFTER  ALL  PIVOT  OPERATIONS  HAVE  BEEN  PERFORMED 
C 

600  CO  710  1  =  1, N 
ll=M+l-I 

K=INDEX( II 1/4096 
ICOLUM=INDEX< U )-4096*K 
IF  (K.EQ. ICCLUM)  GO  TO  710 
CO  705  J=1,N 
SWAP=A( J,K) 
A( J,K)=A( J, ICOLUM) 
A( J,  1C0LUM)=SWAP 
705  CONTINUE 

710  CONTINUE 
C 

711  CONTINUE 
RETURN 
END 


ANGLE 

FUNCTION  ANGLE1RC.RS) 

IMPLICIT  REAL*8(A-U,C-Z) 

K=DSCRT(RC**2+RS**2) 

ANGLE=DARCOS(RC/R) 

IFIRC/R.GT. .900)  ANGLE=  DS IGM ( ANGLE ,RS) 

IFIRC/R.LT.  .900  .AND.  RS.LT.  O.DO)  ANGLE=6.283 1 853-ANGLE 

RETURN 

END 
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ABSTRACT 

The  temperature  wave  method  of  determining  the  thermal  diffusivity 
and  the  thermal  contact  conductance  was  investigated.   The  theoretical 
development  of  the  plane  temperature  wave  behavior  at  the  contact  of 
two  solids  is  given  for  an  interface  perpendicular  to  the  direction  of 
wave  travel.   The  relative  amplitude  and  the  phase  angle  of  the  reflected 
and  the  transmitted  wave  in  relation  to  the  approaching  wave  were 
calculated  in  relation  to  the  dependence  of  the  thermal  properties  of 
the  two  contacting  solids  and  the  thermal  conductance  of  the  contact. 

The  experimental  apparatus  is  described.   The  temperature  waves 
were  generated  by  alternately  heating  and  cooling  one  end  of  the  specimen 
(specimen  pair)  while  the  other  end  was  kept  at  a  constant  temperature. 
The  temperatures,  represented  by  the  thermocouple  voltages,  were  recorded 
on  a  magnetic  tape  recorder.   The  records  were  analyzed  using  a  hybrid 
analog-digital  correlation  computer.   The  theory  of  the  correlation 
function  is  outlined  and  a  computer  program  is  described  which  performs  a 
least  squares  fitting  of  a  set  of  data  with  an  analytical  model  containing 
one  nonlinear  parameter.   The  program  is  particularly  adapted  for  the 
analysis  of  the  correlation  function  containing  a  prominent  harmonic  and 
compensation  for  some  imperfections  of  the  electronic  equipment  used  are 
included. 

Results  are  given  of  the  measurements  of  the  thermal  diffusivity  of 
the  aluminum  alloy  2024-T4  and  Armco  iron  and  of  the  thermal  contact 
conductance  of  the  aluminum-aluminum  and  aluminum-iron  contact  in  air 


at  different  applied  pressures.  The  agreement  with  the  results  from 
literature  is  very  good  for  the  thermal  diffusivity  measurements  and 
good  for -the  thermal  contact  conductance  measurements. 

The  temperature  wave  method  was  found  adequate  for  measuring  the 
thermal  contact  conductance  and  suggestions  are  made  for  some  further 
interesting  investigations  using  this  technique. 


